Revista Alma Mater
Revista Alma Mater
e-ISSN: 3060-9917
2026 | Revista Alma Mater, 1(1), 1–37
www.revistaalmamater.org

Beam Models and Quantum Systems with Linear Potentials Based on Airy Functions

Modelos de vigas y sistemas cuánticos con potenciales lineales basados en funciones de Airy

Juan Toribio Milane
Institute of Mathematics and Institute of Physics, Faculty of Sciences, Universidad Autónoma de Santo Domingo (UASD), Dominican Republic
José Angel Gómez Hernández*
Institute of Mathematics, Faculty of Sciences, Universidad Autónoma de Santo Domingo (UASD), Dominican Republic
Francis Leandro Álvarez Paulino
Institute of Mathematics, Faculty of Sciences, Universidad Autónoma de Santo Domingo (UASD), Dominican Republic
Manuel Leonardo Reyes Cordero
Department of Teaching, Universidad Nacional Evangélica (UNEV), Dominican Republic

*Corresponding author: José Angel Gómez Hernández; jgomez96@uasd.edu.do

Received: 20/09/2025  |  Accepted: 28/03/2026  |  Published: 14/04/2026
DOI: https://doi.org/10.61780/84rzpz60
Cómo citar: Toribio Milane, J., Gómez Hernández, J. A., Álvarez Paulino, F. L. y Reyes Cordero, M. L. (2026). Beam Models and Quantum Systems with Linear Potentials Based on Airy Functions. Revista Alma Mater, 1(1), 1–37. https://doi.org/10.61780/84rzpz60

Beam Models and Quantum Systems with Linear Potentials Based on Airy Functions.
Modelos de vigas y sistemas cuánticos con potenciales lineales basados en funciones de Airy

Institute of Mathematics and Institute of Physics, Faculty of Sciences, Universidad Autónoma de Santo Domingo (UASD), Dominican Republic
ORCID: 0000-0002-0782-1827
Institute of Mathematics, Faculty of Sciences,Universidad Autónoma de Santo Domingo (UASD), Dominican Republic
ORCID: 0009-0009-8182-5840
Institute of Mathematics, Faculty of Sciences, Universidad Autónoma de Santo Domingo (UASD), Dominican Republic
ORCID: 0009-0001-3712-9516
Department of Teaching, Universidad Nacional Evangélica (UNEV), Dominican Republic
ORCID:0009-0000-5972-962X

Received: 20/09/2025 Accepted: 28/03/2026 Published: 14/04/2026
Corresponding author: José Angel Gómez Hernández; jgomez96@uasd.edu.do

Abstract. This paper presents a unified analytical treatment of Airy functions across three canonical models with linear-potential structure: the Euler–Bernoulli beam under self-weight, the quantum bouncer, and a charged particle in a uniform electric field. Starting from the Airy equation, we derive the corresponding series, Bessel-type, and contour-integral representations, and show that admissibility, spectral conditions, and characteristic scales are governed by a common Airy-function mechanism. In the beam problem, a Bessel-type condition determines the critical buckling lengths, whereas in the quantum models quantization is governed by the zeros of \(\operatorname{Ai}\). These results reveal a common analytical framework connecting special-function theory with representative problems in continuum and quantum mechanics.

Keywords: Airy functions, Bessel functions, Euler–Bernoulli beam, quantum bouncer, uniform electric field.

Resumen. Este trabajo presenta un tratamiento analítico unificado de las funciones de Airy en tres modelos canónicos con estructura de potencial lineal: la viga de Euler–Bernoulli sometida a su propio peso, el rebote cuántico y una partícula cargada en un campo eléctrico uniforme. Partiendo de la ecuación de Airy, se derivan representaciones en series, de tipo Bessel y mediante integrales de contorno, y se demuestra que la admisibilidad, las condiciones espectrales y las escalas características están regidas por un mecanismo común basado en las funciones de Airy. En el problema de la viga, una condición de tipo Bessel determina las longitudes críticas de pandeo, mientras que en los modelos cuánticos la cuantización está gobernada por los ceros de \(\operatorname{Ai}\). Estos resultados evidencian un marco analítico común que vincula la teoría de funciones especiales con problemas representativos de la mecánica del continuo y la mecánica cuántica.

Palabras clave: funciones de Airy, funciones de Bessel, viga de Euler–Bernoulli, rebotador cuántico, campo eléctrico uniforme.

10pt 10pt . 0.5em

Introduction

Airy functions play a central role in mathematical physics, as they provide an exact description of the transition between oscillatory and exponential regimes. This property makes them fundamental in the analysis of boundary and turning-point phenomena in quantum mechanics, optics, wave propagation, and cosmology (, , ).

Their importance is further underscored by several active research directions. In models of wave-like dark matter, interference patterns near gravitational caustics are accurately described by the Airy function, and the fringe spacing observed in Schrödinger–Poisson simulations can be explained through the zeros of \(\mathrm{Ai}\), leading to astrophysically relevant predictions . In few-body quantum systems, the zeros of \(\mathrm{Ai}\) and \(\mathrm{Ai}'\) determine leading corrections to the small-mass energy spectrum, highlighting the role of Airy functions in spectral quantization and in transitions between different quantum statistics . Similarly, modified Airy-based semiclassical methods significantly improve standard WKB approximations near turning points, particularly in the analysis of double-well potentials and potential barriers .

From the standpoint of special-function theory, the Airy family continues to motivate new analytical developments. Explicit expressions for products and derivatives of \(\mathrm{Ai}\) and \(\mathrm{Bi}\), together with their connections to generalized hypergeometric functions, reinforce their structural relevance in modern analysis and computational mathematics . Beyond this theoretical framework, Airy functions also play a significant role in optics and wave engineering. In particular, Airy-type beams have become an important tool for manipulating structured waves and quantum states of light due to their non-diffracting and self-accelerating properties . This line of research is complemented by the experimental realization of electron Airy beams in microscopy , as well as by extensions involving Airy–Gaussian beams in paraxial systems.

The mathematical theory itself is classical and well documented. Comprehensive treatments of Airy functions and their role within special-function theory can be found in standard references such as the NIST Handbook of Mathematical Functions , the monograph by Vallée and Soares , and the tables of Gradshteyn and Ryzhik . Their significance in quantum mechanics is equally well established: Landau and Lifshitz and Griffiths and Schroeter show that the Schrödinger equation with a linear potential naturally reduces to the Airy equation, while the associated turning-point analysis is closely linked to the WKB connection formalism . These references make clear that the Airy function is not an isolated construct within special-function theory, but a recurring analytical mechanism in canonical physical models. This background, however, also reveals a persistent fragmentation. The analytic construction of Airy functions, their connection with Bessel equations, their contour-integral representations, and their role in concrete continuum and quantum models are typically developed in separate contexts, with differing emphases and without an explicit derivational synthesis. Consequently, the common analytical mechanism underlying these problems often remains implicit.

The aim of this paper is to make this mechanism explicit. Rather than presenting a general survey of Airy-function theory, we develop a focused analytical treatment of the representations and reductions required to study three canonical models governed by linear potentials or equivalent linearized forces: the Euler–Bernoulli beam under self-weight, the quantum bouncer, and a charged particle in a uniform electric field. In this framework, Airy functions are not merely auxiliary solutions; they determine the admissible branches, the asymptotic selection mechanism, the spectral quantization rules, and the characteristic physical scales of the models under consideration.

This perspective is particularly natural in quantum systems with linear potentials. A representative example arises in semiconductor physics, where band bending under an approximately uniform electric field leads directly to the Airy equation and . Similar reductions occur in gravitational and electrostatic one-dimensional models, where the zeros of \(\operatorname{Ai}\) encode the discrete spectrum. What is common to these settings is not only the same differential equation, but also the same special-function mechanism governing boundary admissibility and physically meaningful quantization.

This viewpoint also clarifies the scope of the present work relative to standard references. In , Airy functions arise in the analysis of triangular wells in semiconductor physics; in , they appear in the classical reduction of Schrödinger’s equation with a linear potential; and in the quantum bouncer is treated as a specific spectral model. What the present paper adds is an explicit analytical synthesis that places the series construction, Bessel reductions, contour-integral formulation, and zero-driven spectral conditions within a unified framework encompassing both continuum and quantum models.

Our contribution is therefore one of rigorous synthesis and unified exposition at the level of analytical structure. More precisely, the manuscript does not aim to survey the full Airy-function literature, but to assemble, within a single derivational framework, the ingredients essential for the treatment of three canonical models with linear-potential structure. We first derive the Airy solutions via power series and fix the classical normalizations at \(z=0\) through the Gamma function. We then establish their connection with ordinary and modified Bessel functions by means of changes of variables that separate the exponential (\(x>0\)) and oscillatory (\(x<0\)) regimes. Next, we develop the contour-integral representation \(y(z)=\int_C f(t)e^{zt}\,dt\), identify the decay sectors \(\Re(t^3)>0\), determine the canonical contours \(C_1,C_2,C_3\), and justify the emergence of the independent solutions \(\operatorname{Ai}\) and \(\operatorname{Bi}\). These tools are subsequently applied to the three models above, where we derive the corresponding critical or quantization conditions and show how the relevant scales arise directly from the Airy structure itself.

The paper is organized as follows. Section 2 develops the analytical framework for Airy functions, including their series construction, Bessel reductions, contour representations, and the distribution of their zeros. Section 3 applies this framework to the Euler–Bernoulli beam under self-weight and to two one-dimensional quantum models, deriving the associated spectral conditions and characteristic scales. Section 4 discusses the analytical contribution of the paper in relation to the existing literature, and Section 5 presents the conclusions.

Preliminaries

We consider the Airy differential equation \[\label{eq:Airy} y''(x) - x\,y(x) = 0 ,\] whose coefficients are entire functions on \(\mathbb{C}\). By the classical Cauchy–Kovalevskaya theorem for linear ordinary differential equations with analytic coefficients (see, e.g., , Lec. 4), every solution of [eq:Airy] admits a convergent power series expansion around any prescribed point \(x_0\in\mathbb{R}\). Without loss of generality, we set \(x_0=0\) and employ the series ansatz \[\label{eq:series-ansatz} y(x) = \sum_{m=0}^{\infty} c_m x^m .\]

Differentiating term by term, we obtain \[y'(x)=\sum_{m=0}^{\infty} (m+1)c_{m+1}x^m, \qquad y''(x)=\sum_{m=0}^{\infty} (m+1)(m+2)c_{m+2}x^m,\] whereas multiplication by \(-x\) yields \[-\,x\,y(x)=-\sum_{m=1}^{\infty} c_{m-1} x^{m}.\]

Substituting these series into [eq:Airy] and equating coefficients of like powers of \(x\), we obtain the recurrence relations \[\sum_{m=0}^{\infty} \bigl[(m+1)(m+2)c_{m+2}\bigr]x^m -\sum_{m=1}^{\infty} c_{m-1}x^m = 0.\]

By isolating the case \(m=0\) and comparing coefficients, one derives \[\label{eq:airy-recurrence} 2c_2=0,\qquad (m+1)(m+2)c_{m+2}=c_{m-1},\quad m\ge1.\] Thus \(c_2=0\), and the recurrence [eq:airy-recurrence] links \(c_{m+2}\) with \(c_{m-1}\). Consequently, the coefficients separate into three disjoint arithmetic progressions according to their indices modulo \(3\). In particular, the congruence class \(m\equiv 2 \pmod{3}\) vanishes identically, i.e., \[\label{eq:c3k+2zero} c_{3k+2}=0,\qquad k\ge0.\]

From [eq:airy-recurrence] one deduces explicit formulas for the subsequences. For even multiples of three, \[\label{eq:c3k} c_{3k}=\frac{1}{(3k)!}\prod_{i=1}^{k}(3i-2)\,c_0,\qquad k\ge1,\] while for odd multiples of three plus one, \[\label{eq:c3k+1} c_{3k+1}=\frac{1}{(3k+1)!}\prod_{i=1}^{k}(3i-1)\,c_1,\qquad k\ge1,\] with \(c_0\) and \(c_1\) arbitrary. Therefore the general solution decomposes as \[\label{eq:airy-basis} y(x)=c_0\,y_1(x)+c_1\,y_2(x),\] where \[y_1(x)=1+\sum_{k=1}^{\infty}\frac{1}{(3k)!}\prod_{i=1}^{k}(3i-2)\,x^{3k},\quad y_2(x)=x+\sum_{k=1}^{\infty}\frac{1}{(3k+1)!}\prod_{i=1}^{k}(3i-1)\,x^{3k+1}.\] At \(x=0\) these series satisfy \(y_1(0)=1,\,y_1'(0)=0\) and \(y_2(0)=0,\,y_2'(0)=1\), so that \(W(y_1,y_2)(0)=1\neq0\). Thus \(y_1,y_2\) are linearly independent and provide a fundamental system of solutions, see .

The convergence of both series is established by the ratio test. For example, if \(T_k(x)\) denotes the \(k\)–th summand of \(y_1\), then \(|T_{k+1}/T_k|\to0\) as \(k\to\infty\), and similarly for the summands \(S_k(x)\) of \(y_2\). Consequently, each series has infinite radius of convergence. A uniform bound on compact subsets (cf. Remark 1, see also and ) shows that all derivatives also converge uniformly, so \(y_1\) and \(y_2\) are entire functions.

Remark 1 (2). For any \(R>0\), there exists \(K(R)\) such that the ratios of consecutive terms of \(y_1\) and \(y_2\) are bounded by \(1/2\) whenever \(k \ge K(R)\) and \(|x| \le R\). By the Weierstrass \(M\)–test, both series converge uniformly on compact disks, which justifies differentiation term by term (see, for instance, (, Ch. 8)).

It follows that the classical Airy functions can be expressed as the normalized combinations \[\label{eq:AiBi-defi} \operatorname{Ai}(x)=\frac{1}{3^{2/3}\Gamma(\tfrac{2}{3})}\,y_1(x)-\frac{1}{3^{1/3}\Gamma(\tfrac{1}{3})}\,y_2(x),\,\, \operatorname{Bi}(x)=\frac{1}{3^{1/6}\Gamma(\tfrac{2}{3})}\,y_1(x)+\frac{3^{1/6}}{\Gamma(\tfrac{1}{3})}\,y_2(x),\] where \(\Gamma\) denotes Euler’s gamma function. These choices enforce the initial conditions \[\label{eq:AiBi-initiali} \operatorname{Ai}(0)=\frac{1}{3^{2/3}\Gamma(\tfrac{2}{3})},\quad \operatorname{Ai}'(0)=-\frac{1}{3^{1/3}\Gamma(\tfrac{1}{3})},\quad \operatorname{Bi}(0)=\frac{1}{3^{1/6}\Gamma(\tfrac{2}{3})},\quad \operatorname{Bi}'(0)=\frac{3^{1/6}}{\Gamma(\tfrac{1}{3})},\] so that \(\operatorname{Ai},\operatorname{Bi}\) form the canonical basis of solutions of [eq:Airy]. As \(y_1,y_2\) are entire, the combinations [eq:AiBi-defi] show that \(\operatorname{Ai}\) and \(\operatorname{Bi}\) are also entire and analytic throughout the complex plane.

Bessel equations and their connection with the Airy equation

The ordinary Bessel functions \(J_\nu\) and \(Y_\nu\) satisfy \[\label{eq:bessel-clasica} x^2 u''(x) + x u'(x) + \bigl(x^2 - \nu^2\bigr) u(x) = 0,\] so that the general solution is \(u(x) = A J_{\nu}(x) + B Y_{\nu}(x)\), where \[J_{\nu}(x) = \sum_{k=0}^{\infty} \frac{(-1)^{k}}{k!\,\Gamma(k+\nu+1)} \left(\frac{x}{2}\right)^{2k+\nu}, \quad Y_{\nu}(x) = \frac{ J_{\nu}(x)\cos(\nu \pi) - J_{-\nu}(x) }{ \sin(\nu \pi) }, \quad \forall \nu \notin \mathbb{Z}.\] In the case of integer orders \(n \in \mathbb{Z}\), the function is defined as the limit of the family \(\{ Y_{\nu}(x)\}_{\nu \notin \mathbb{Z}}\) as \(\nu \to n\): \[Y_{n}(x) = \lim_{\nu \to n} Y_{\nu}(x), \qquad n \in \mathbb{Z}.\]

Similarly, the modified Bessel functions \(I_\nu\) and \(K_\nu\) satisfy \[\label{eq:bessel-modificada} x^2 u''(x) + x u'(x) - \bigl(x^2 + \nu^2\bigr) u(x) = 0,\] with general solution \(u(x) = A I_{\nu}(x) + B K_{\nu}(x)\), where \[I_{\nu}(x) = \sum_{k=0}^{\infty} \frac{1}{k!\,\Gamma(k+\nu+1)} \left( \frac{x}{2} \right)^{2k+\nu}, \qquad K_{\nu}(x) = \frac{\pi}{2}\, \frac{ I_{-\nu}(x) - I_{\nu}(x) }{ \sin(\nu \pi) }, \quad \forall \nu \notin \mathbb{Z}.\] For integer order \(n \in \mathbb{Z}\), the function is defined as the following limit: \[K_{n}(x) = \lim_{\nu \to n} K_{\nu}(x), \qquad n \in \mathbb{Z}.\]

See , .

Proposition 1 (Transformation for \(x>0\)). Let \(x>0\) and set \[\xi=\tfrac{2}{3}x^{3/2},\qquad y(x)=x^{1/2}u(\xi).\] Then \(y\) satisfies [eq:Airy] if and only if \(u\) satisfies [eq:bessel-modificada] with order \(\nu=\tfrac13\), namely \[\label{eq:bessel-u-positiva} u''(\xi)+\tfrac{1}{\xi}u'(\xi)-\Bigl(1+\tfrac{1}{9\xi^2}\Bigr)u(\xi)=0.\] Consequently, the general solution of [eq:Airy] for \(x>0\) is \[\label{eq:sol-general-xpos} y(x)=\sqrt{x}\,\Bigl\{C_1 I_{1/3}(\xi)+C_2 I_{-1/3}(\xi)\Bigr\} =\sqrt{x}\,\Bigl\{A K_{1/3}(\xi)+B I_{1/3}(\xi)\Bigr\}.\]

Proof. With \(y=x^{1/2}u(\xi)\) and \(\xi'=\sqrt{x}\), we compute \[y'=\tfrac12 x^{-1/2}u+x\,u',\qquad y''=-\tfrac14 x^{-3/2}u+\tfrac32 u'+x^{3/2}u'',\] so that \(y''-xy=0\) reduces to \[\xi^2 u''+\xi u'-(\xi^2+\tfrac19)u=0,\] which is precisely [eq:bessel-modificada] with \(\nu=\tfrac13\). ◻

Proposition 2 (Transformation for \(x<0\)). Let \(x=-s<0\) with \(s>0\) and set \[\eta=\tfrac{2}{3}s^{3/2},\qquad y(-s)=s^{1/2}v(\eta).\] Then \(y\) satisfies [eq:Airy] if and only if \(v\) satisfies [eq:bessel-clasica] with \(\nu=\tfrac13\), that is, \[\label{eq:bessel-v-negativa} v''(\eta)+\tfrac{1}{\eta}v'(\eta)+\Bigl(1-\tfrac{1}{9\eta^2}\Bigr)v(\eta)=0.\] Therefore, the general solution for \(x<0\) is \[\label{eq:sol-general-xneg} y(x)=\sqrt{|x|}\,\Bigl\{D_1 J_{1/3}(\eta)+D_2 J_{-1/3}(\eta)\Bigr\} =\sqrt{|x|}\,\Bigl\{E_1 J_{1/3}(\eta)+E_2 Y_{1/3}(\eta)\Bigr\}.\]

Using \(\xi=\tfrac{2}{3}x^{3/2}\) for \(x>0\) and \(\eta=\tfrac{2}{3}|x|^{3/2}\) for \(x<0\), together with the values at \(x=0\), \[\operatorname{Ai}(0)=\tfrac{1}{3^{2/3}\Gamma(\tfrac23)},\quad \operatorname{Ai}'(0)=-\tfrac{1}{3^{1/3}\Gamma(\tfrac13)},\qquad \operatorname{Bi}(0)=\tfrac{1}{3^{1/6}\Gamma(\tfrac23)},\quad \operatorname{Bi}'(0)=\tfrac{3^{1/6}}{\Gamma(\tfrac13)},\] and using the expansions of Bessel functions near \(z=0\) , one derives the connection formulas.

For positive arguments, the Airy functions can be expressed in terms of modified Bessel functions: \[\begin{aligned} \operatorname{Ai}(x)&=\tfrac{1}{3}\sqrt{x}\,\bigl[I_{-1/3}(\xi)-I_{1/3}(\xi)\bigr] =\tfrac{1}{\pi}\sqrt{\tfrac{x}{3}}\,K_{1/3}(\xi), \label{secondequality} \\[2mm] \operatorname{Bi}(x)&=\sqrt{\tfrac{x}{3}}\,\bigl[I_{-1/3}(\xi)+I_{1/3}(\xi)\bigr]. \end{aligned}\] Their asymptotics as \(x\to+\infty\) reflect exponential decay and growth: \[\label{eq:Ai-Bi-asymp-plus} \operatorname{Ai}(x)\sim \tfrac{1}{2\sqrt{\pi}}x^{-1/4}e^{-\tfrac{2}{3}x^{3/2}},\qquad \operatorname{Bi}(x)\sim \tfrac{1}{\sqrt{\pi}}x^{-1/4}e^{+\tfrac{2}{3}x^{3/2}}.\]

For negative arguments, the reduction involves ordinary Bessel functions, giving \[\begin{aligned} \operatorname{Ai}(x)&=\tfrac{1}{3}\sqrt{|x|}\,\bigl[J_{-1/3}(\eta)+J_{1/3}(\eta)\bigr], \label{eq:Ai-negativa-bessel}\\[2mm] \operatorname{Bi}(x)&=\tfrac{1}{\sqrt{3}}\sqrt{|x|}\,\bigl[J_{-1/3}(\eta)-J_{1/3}(\eta)\bigr]. \label{eq:Bi-negativa-bessel} \end{aligned}\] The oscillatory character is explicit in the asymptotics as \(x\to -\infty\): \[\label{eq:Ai-Bi-asymp-minus} \operatorname{Ai}(x)\sim \tfrac{1}{\sqrt{\pi}}|x|^{-1/4} \sin\!\Bigl(\tfrac{2}{3}|x|^{3/2}+\tfrac{\pi}{4}\Bigr),\qquad \operatorname{Bi}(x)\sim \tfrac{1}{\sqrt{\pi}}|x|^{-1/4} \cos\!\Bigl(\tfrac{2}{3}|x|^{3/2}+\tfrac{\pi}{4}\Bigr).\]

Finally, the second equality in [secondequality] follows from the identity \[\label{eq:K-identity} K_\nu(z)=\frac{\pi}{2\sin(\nu\pi)}\bigl(I_{-\nu}(z)-I_{\nu}(z)\bigr),\] valid for \(\nu=\tfrac13\) .

Integral representations of \(\operatorname{Ai}\) and \(\operatorname{Bi}\)

Consider the Airy equation in its canonical form: \[\label{eq:airy} y''(z)=z\,y(z), \qquad z\in\mathbb{C}.\]

We look for solutions as contour integrals in the \(t\)–plane (see ): \[\label{eq:ansatz} y(z)=\int_{C} f(t)\,e^{z t}\,dt,\] where \(C\) is a piecewise oriented contour in the \(t\)–plane and \(f\) is a function to be determined. We say that \(C\) is admissible if its endpoints \(a,b\) lie at infinity along rays, and it is smooth or at least piecewise differentiable, \[t=r e^{i\theta},\quad r\to+\infty,\] contained in angular sectors where \(\Re(t^3)>0\); equivalently, \(\cos(3\theta)>0\), that is, \[\label{eq:sectores} 2\pi n-\tfrac{\pi}{2}<3\arg t<2\pi n+\tfrac{\pi}{2},\qquad n\in\mathbb{Z},\] which determines exactly three open sectors of width \(\pi/3\) centered at \(\arg t=0,\,2\pi/3,\,4\pi/3\). We also assume that \(C\) has finite length inside compact regions and that [eq:ansatz] converges absolutely.
Let \(t=r e^{i\theta}\). For large \(r\), \[\Re\!\Big(-\tfrac{1}{3}t^3+z t\Big)=-\tfrac{r^3}{3}\cos(3\theta)+r|z|\cos(\theta-\arg z).\] There exist \(c>0\) and \(R_0\) such that \(\Re(-t^3/3+zt)\le -c\,r^3\) for \(r\ge R_0\) along the ends of \(C\). In particular, for any integer \(m\ge 0\), \[\int_{C}\! |t|^m e^{\Re(-t^3/3+z t)}|dt|<\infty,\] and by the dominated convergence theorem, \(y\) is entire in \(z\) and may be differentiated twice under the integral sign: \[\label{eq:derivs} y'(z)=\int_C t\,f(t)\,e^{z t}\,dt,\qquad y''(z)=\int_C t^2 f(t)\,e^{z t}\,dt.\]

From [eq:ansatz], \[z\,y(z)=\int_C z f(t)\,e^{z t}\,dt =\int_C \Big(\tfrac{d}{dt}\big(f(t)e^{z t}\big)-f'(t) e^{z t}\Big)dt =\Big[f(t)e^{z t}\Big]_{a}^{b}-\int_C f'(t)e^{z t}\,dt.\] If \(C\) is admissible and \(f(t)=\mathcal{O}(e^{-\Re(t^3)/3})\), then \(\big[f(t)e^{z t}\big]_{a}^{b}=0\), since \(|f(t)e^{z t}|\to 0\) at the contour ends. Using [eq:ansatz], the Airy equation [eq:airy] reduces to \[\int_C \big(t^2 f(t)+f'(t)\big)e^{z t}\,dt=0 \quad \text{for all } z\in\mathbb{C}.\] Since the integrand is entire in \(t\) and \(e^{z t}\) never vanishes, we obtain the differential equation for \(f\): \[\label{eq:ode_f} f'(t)+t^2 f(t)=0.\]

The equation [eq:ode_f] is first order linear, and its general solution is \[\label{eq:f_sol} f(t)=A\,\exp\!\Big(-\tfrac{t^3}{3}\Big),\qquad A\in\mathbb{C}.\] Note that \(t\mapsto t^3\) is an entire polynomial, so [eq:f_sol] does not require branch choices.

Thus, for any admissible contour \(C\), the function \[\label{eq:airy_integral_form} y(z)=A\int_C \exp\!\Big(-\tfrac{t^3}{3}+z t\Big)\,dt\] is well defined, entire in \(z\), and satisfies [eq:airy]. The constant \(A\) and the choice of the pair of sectors where the contour ends lie determine a family of solutions; specific choices of \(C\) yield the classical Airy functions.
We now analyze the decay sectors of the integrand. For \(t=re^{i\theta}\), the modulus in [eq:airy_integral_form] is \[\big|e^{-t^{3}/3+zt}\big|=\exp\!\Big(-\tfrac{r^{3}}{3}\cos(3\theta)+r|z|\cos(\theta-\arg z)\Big).\] As \(r\to\infty\), the cubic term dominates; thus exponential decay occurs along rays with \[\Re(t^{3})>0\quad\Longleftrightarrow\quad \cos(3\theta)>0 \quad\Longleftrightarrow\quad 2\pi n-\tfrac{\pi}{2}<3\theta<2\pi n+\tfrac{\pi}{2},\quad n\in\mathbb Z.\] This determines exactly three open sectors of aperture \(\pi/3\): \[S_k=\bigl\{t\neq 0:\ |\arg t-\tfrac{2\pi k}{3}|<\tfrac{\pi}{6}\bigr\},\qquad k=0,1,2,\] centered at the angles \(0,\,2\pi/3,\,4\pi/3\), respectively. A pair of endpoints is admissible if each lies at infinity in distinct sectors among \(\{S_0,S_1,S_2\}\).

In this case, by the previous estimate, there exists \(c>0\) such that \[\Big|\exp\!\big(-\tfrac{t^{3}}{3}+zt\big)\Big|\le e^{-c r^{3}}\,e^{O(r)}\qquad(r\to\infty \text{ along the contour ends}),\] so that [eq:airy_integral_form] converges absolutely; moreover, the boundary term in the integration by parts vanishes.
If both ends lie in the same sector, the integral is zero. If both belong to the same \(S_k\), joining them by a large circular arc \[\Gamma_R=\{R e^{i\phi}:\ \phi\in[\theta_1,\theta_2]\subset S_k\},\] and considering the closed loop \(C_R:=C\cup\Gamma_R\), since the integrand is entire we have \(\int_{C_R}e^{-t^3/3+zt}\,dt=0\). Furthermore, \[\Big|\int_{\Gamma_R}\!e^{-t^3/3+zt}\,dt\Big|\le \int_{\Gamma_R}\!\exp\!\big(-\tfrac{R^3}{3}\cos(3\phi)+O(R)\big)\,|dt| \le (\theta_2-\theta_1)R\,e^{-cR^3+O(R)}\xrightarrow[R\to\infty]{}0,\] since \(\cos(3\phi)\ge\delta>0\) inside the sector. Hence, in the limit \(R\to\infty\), one has \(\int_C e^{-t^3/3+zt}\,dt=0\).
As the ends must belong to distinct sectors, there are exactly three canonical connections: \[C_1:\ S_2\longrightarrow S_1, \quad C_2:\ S_0\longrightarrow S_1, \quad C_3:\ S_0\longrightarrow S_2,\] that is, joining the central rays \(4\pi/3\to 2\pi/3\), \(0\to 2\pi/3\), and \(0\to 4\pi/3\), respectively. The orientation (from the first sector listed to the second) fixes the sign of each integral (see Fig. 1).

Decay sectors \(\Re(t^3)>0\) in the \(t\)–plane and the canonical contours \(C_1,C_2,C_3\) used in [eq:Ai-def][eq:Bi-def].

Let \[I(C)=\int_C e^{-t^3/3+zt}\,dt.\] For large \(R\), consider the polygonal contour \(L_R\) joining \(R e^{i\cdot 0}\), \(R e^{i\cdot 2\pi/3}\) and \(R e^{i\cdot 4\pi/3}\) by internal arcs contained within the decay sectors \(S_k\). Then the boundary of the region delimited by \(L_R\) can be traversed, in the limit \(R\to\infty\), as \[C_2 + C_3 - C_1,\] and since the integrand is entire, \(\lim_{R\to\infty}\int_{L_R} e^{-t^3/3+zt}\,dt=0\). Thus \[I(C_2)+I(C_3)-I(C_1)=0,\] which shows that only two integrals are independent. This justifies that it suffices to choose two convenient linear combinations to define the special Airy functions.

Let \(\Omega=e^{2\pi i/3}\). The rotations \(t\mapsto\Omega t\) permute the sectors \(S_k\) and, with appropriate orientation, map \(C_1,C_2,C_3\) into each other. Consequently, if \(I_k(z):=I(C_k)\), then \[I_1(z)+\Omega\, I_1(\Omega z)+\Omega^2 I_1(\Omega^2 z)=0,\] and linear combinations of the \(I_k\) with rotated arguments allow all solutions to be expressed in terms of two integrals.
As the Airy equation is second order, two linearly independent solutions are required. The classical convention, used in the standard literature , is to define: \[\label{eq:Ai-def} \operatorname{Ai}(z)=\frac{1}{2\pi i}\int_{C_1}\exp\!\Big(-\tfrac{t^3}{3}+zt\Big)\,dt,\] \[\label{eq:Bi-def} \operatorname{Bi}(z)=\frac{1}{2\pi}\int_{C_2}\exp\!\Big(-\tfrac{t^3}{3}+zt\Big)\,dt -\frac{1}{2\pi}\int_{C_3}\exp\!\Big(-\tfrac{t^3}{3}+zt\Big)\,dt.\]

The factors \(1/(2\pi i)\) and \(1/(2\pi)\) are chosen so that:

These definitions are justified as follows:

  1. The choice of \(C_1\) in [eq:Ai-def] is natural, since deforming this contour towards the imaginary axis yields the real representation \[\label{airyreal} \operatorname{Ai}(z)=\frac{1}{\pi}\int_{0}^{\infty}\cos\!\Big(\tfrac{s^3}{3}+z s\Big)\,ds,\]

    which shows that \(\operatorname{Ai}(x)\in\mathbb{R}\) for \(x\in\mathbb{R}\) and that \(\operatorname{Ai}(x)\to 0\) as \(x\to+\infty\).

  2. The definition of \(\operatorname{Bi}(z)\) in [eq:Bi-def] symmetrically combines integrals over \(C_2\) and \(C_3\), which ensures that \(\operatorname{Bi}(x)\) is real for \(x\in\mathbb{R}\), and that for \(x\to-\infty\) its asymptotic behavior is \[\label{asintoticbi} \operatorname{Bi}(x)\sim \frac{1}{\sqrt{\pi}|x|^{1/4}}\cos\!\Big(\tfrac{2}{3}|x|^{3/2}+\tfrac{\pi}{4}\Big),\] while \[\label{asintoticai} \operatorname{Ai}(x)\sim \frac{1}{\sqrt{\pi}|x|^{1/4}}\sin\!\Big(\tfrac{2}{3}|x|^{3/2}+\tfrac{\pi}{4}\Big),\] that is, \(\operatorname{Bi}\) is in quadrature with \(\operatorname{Ai}\) (see, e.g., ).

Let \(\Omega=e^{2\pi i/3}\). The invariance of the Airy equation under \(z\mapsto\Omega z\) and the action of \(\Omega\) on the contours produce the identities: \[\label{eq:rot-id} \operatorname{Ai}(z)+\Omega\,\operatorname{Ai}(\Omega z)+\Omega^{2}\,\operatorname{Ai}(\Omega^{2}z)=0,\] \[\label{eq:Bi-rot} \operatorname{Bi}(z)=i\big(\Omega^{2}\operatorname{Ai}(\Omega^{2}z)-\Omega\,\operatorname{Ai}(\Omega z)\big) =e^{-\pi i/6}\operatorname{Ai}(ze^{-2\pi i/3})+e^{\pi i/6}\operatorname{Ai}(ze^{2\pi i/3}).\] These relations allow \(\operatorname{Bi}\) to be expressed in terms of \(\operatorname{Ai}\) evaluated at rotated arguments.
Thus, the contour integrals [eq:Ai-def][eq:Bi-def] canonically define the two Airy functions that form a fundamental basis of solutions of [eq:airy].

Finally, consider again the contour definition in [eq:Ai-def], where \(C_1\) connects the decay sectors centered at \(4\pi/3\) and \(2\pi/3\). Our aim is to deform \(C_1\) to (almost) the imaginary axis and deduce the real form from [airyreal] (see Fig. 2)

Deformation of the canonical contour \(C_1\) toward the vertical line \(\Re t=-\varepsilon\) inside the light-green decay sectors \(\Re(t^3)>0\). The inset shows the before/after scheme. Notation follows the text: \(S_k\), \(C_1\), \(\Gamma_{4\pi/3}(R)\), \(\Gamma_{2\pi/3}(R)\), and \(V_\varepsilon\).

Step 1: Deformation of \(C_1\) to a vertical line \(\Re t=-\varepsilon\).

Let \(R>0\) large and \(0<\varepsilon\ll1\). We truncate \(C_1\) to a contour \(C_1(R)\) and deform it (without crossing singularities, since the integrand is entire) into a contour consisting of: (i) an arc \(\Gamma_{4\pi/3}(R)\) inside the decay sector around \(4\pi/3\) descending from \(R e^{4\pi i/3}\) to \(-\varepsilon-iR\); (ii) the vertical segment \[V_\varepsilon(R)=\{-\varepsilon+is:\ s\in[-R,R]\};\] and (iii) an arc \(\Gamma_{2\pi/3}(R)\) inside the decay sector around \(2\pi/3\) ascending from \(-\varepsilon+iR\) to \(R e^{2\pi i/3}\). The orientation goes from \(S_2\) to \(S_1\).

On the arcs \(\Gamma_{\ast}(R)\), since \(\Re(t^3)>0\) inside the decay sectors, we have \[\big|\exp(-t^3/3+zt)\big| =\exp\!\Big(-\tfrac{1}{3}R^3\cos(3\arg t)+O(R)\Big) \le \exp(-cR^3+O(R)),\] for some \(c>0\) independent of \(R\). Hence, \[\int_{\Gamma_{\ast}(R)}\exp(-t^3/3+zt)\,dt \xrightarrow[R\to\infty]{} 0,\] and in the limit \(R\to\infty\) only the contribution of \(V_\varepsilon:=V_\varepsilon(\infty)\) remains. Consequently, \[\label{eq:C1-to-vertical} \operatorname{Ai}(z)=\frac{1}{2\pi i}\int_{V_\varepsilon}\exp\!\Big(-\tfrac{t^3}{3}+zt\Big)\,dt, \qquad (0<\varepsilon\ll1).\]

Step 2: Vertical parametrization and Gaussian damping.

We parametrize \(V_\varepsilon\) by \(t=-\varepsilon+is\), \(s\in\mathbb{R}\), and compute \[\begin{aligned} -\tfrac{t^3}{3}+zt &=-\tfrac{(-\varepsilon+is)^3}{3}+z(-\varepsilon+is) \\ &=\; i\!\left(\tfrac{s^3}{3}+zs\right)\;-\;\varepsilon s^2\;-\;\varepsilon^2 i s\;+\;\tfrac{\varepsilon^3}{3}\;-\;\varepsilon z. \end{aligned}\] Thus, using \(dt=i\,ds\) in [eq:C1-to-vertical], \[\label{eq:eps-form} \operatorname{Ai}(z)=\frac{1}{2\pi}\int_{-\infty}^{\infty} \exp\!\Big(i\!\Big(\tfrac{s^3}{3}+zs\Big)\Big)\, \exp\!\big(-\varepsilon s^2\big)\, \exp\!\big(-\varepsilon^2 i s\big)\, \exp\!\Big(\tfrac{\varepsilon^3}{3}-\varepsilon z\Big)\,ds.\] For \(s\to\pm\infty\), the factor \(e^{-\varepsilon s^2}\) ensures absolute convergence; hence, derivatives with respect to \(z\) may be passed under the integral sign, and dominated convergence applies as \(R\to\infty\).

Step 3: Limit \(\varepsilon\downarrow0\) and odd symmetry of the phase.

Taking \(\varepsilon\downarrow0\) in [eq:eps-form] (as a limit of absolutely convergent integrals) yields the improper oscillatory integral \[\label{eq:Ai-real-symmetric} \operatorname{Ai}(z)=\frac{1}{2\pi}\,\mathrm{p.v.}\!\int_{-\infty}^{\infty} \exp\!\Big(i\!\Big(\tfrac{s^3}{3}+zs\Big)\Big)\,ds,\] where \(\mathrm{p.v.}\) indicates that the limit is taken by symmetric truncation \(\lim_{R\to\infty}\int_{-R}^{R}\exp(i(\tfrac{s^3}{3}+zs))\,ds\).

Note that the phase \(\phi(s)=\tfrac{s^3}{3}+zs\) is an odd function: \(\phi(-s)=-\phi(s)\). Writing \(e^{i\phi(s)}=\cos(\phi(s))+i\sin(\phi(s))\) and truncating on \([-R,R]\), we have that \(\sin\phi\) is odd and \(\cos\phi\) is even; therefore, \[\int_{-R}^{R}\sin\phi(s)\,ds=0,\qquad \int_{-R}^{R}\cos\phi(s)\,ds=2\int_0^{R}\cos\phi(s)\,ds.\] Passing to the limit \(R\to\infty\) gives \[\label{eq:Ai-real-final} \operatorname{Ai}(z)=\frac{1}{\pi}\int_{0}^{\infty} \cos\!\Big(\tfrac{s^3}{3}+zs\Big)\,ds,\] which is precisely the real representation of \(\operatorname{Ai}\) announced in [airyreal].
Finally, to fix the normalization constant in the definition [eq:Ai-def], we evaluate the contour integral of \(\operatorname{Ai}\) at \(z=0\): \[\label{eq:Ai0-int} \operatorname{Ai}(0)=\frac{1}{2\pi i}\int_{C_1}\exp\!\Big(-\tfrac{t^3}{3}\Big)\,dt.\]

Recall that \(C_1\) joins the sectors centered at \(4\pi/3\) and \(2\pi/3\). We may conveniently choose \(C_1\) as the union of two rays: \[\begin{aligned} L_1&:\ t=\rho e^{4\pi i/3},\quad \rho\in[\infty,0],\\ L_2&:\ t=\rho e^{2\pi i/3},\quad \rho\in[0,\infty). \end{aligned}\] That is, \(L_1\) descends from infinity in the direction \(\arg t=4\pi/3\) to the origin, and \(L_2\) ascends from the origin to infinity in the direction \(\arg t=2\pi/3\).
Making the change of variables \(s=\tfrac{1}{3}t^3\), in both segments the mapping \(t\mapsto s=\tfrac{1}{3}t^3\) sends the path to the positive real half–axis (see Fig. 3). Indeed:

Normalization at \(z=0\). The rays \(L_1\) and \(L_2\) defining \(C_1\) are mapped by \(s=\tfrac{1}{3}t^3\) onto the positive real half–axis in the \(s\)–plane, yielding the Gamma–integral used in [eq:Ai0-int]. Notation matches the text: \(L_1\), \(L_2\), and \(\mathbb{R}_+\).

Substituting into [eq:Ai0-int] and collecting contributions: \[\begin{aligned} \operatorname{Ai}(0) &=\frac{1}{2\pi i}\left( \int_{\infty}^{0} e^{-s}\,e^{4\pi i/3}(3s)^{-2/3}ds +\int_{0}^{\infty} e^{-s}\,e^{2\pi i/3}(3s)^{-2/3}ds \right)\\[1ex] &=\frac{3^{-2/3}}{2\pi i}\int_0^{\infty} e^{-s}s^{-2/3} \big(-e^{4\pi i/3}+e^{2\pi i/3}\big)\,ds. \end{aligned}\]

Therefore, \[\operatorname{Ai}(0)=\frac{3^{-2/3}}{2\pi i}\,(i\sqrt{3})\int_0^{\infty}e^{-s}s^{-2/3}\,ds.\]

The resulting integral is the classical definition of the Gamma function: \[\int_0^{\infty} e^{-s}s^{\alpha-1}ds=\Gamma(\alpha),\qquad \Re(\alpha)>0.\] Here \(\alpha-1=-\tfrac{2}{3}\), i.e. \(\alpha=\tfrac{1}{3}\). Thus, \[\int_0^{\infty} e^{-s}s^{-2/3}\,ds=\Gamma\!\Big(\tfrac{1}{3}\Big).\]

Putting everything together: \[\label{eq:Ai0-final} \operatorname{Ai}(0)=\frac{3^{-1/6}}{2\pi}\,\Gamma\!\Big(\tfrac{1}{3}\Big).\]

An analogous argument yields \(\operatorname{Ai}'(0)=-\tfrac{3^{-5/6}}{2\pi}\,\Gamma(\tfrac{2}{3})\), further emphasizing the deep connection between Airy functions and fractional values of the Gamma function.

Zeros of the Airy functions

In the results section we present a quantum model that exploits the zeros of the Airy functions (see Section 3.2). We denote the negative real zeros by \[\mathrm{Ai}(-a_k)=0,\quad \mathrm{Ai}'(-a'_k)=0,\quad \mathrm{Bi}(-b_k)=0,\quad \mathrm{Bi}'(-b'_k)=0,\] with the orderings \[0<a_1<a_2<\cdots,\quad 0<a'_1<a'_2<\cdots,\quad 0<b_1<b_2<\cdots,\quad 0<b'_1<b'_2<\cdots.\]

On the real axis, \(\mathrm{Ai}\), \(\mathrm{Ai}'\), \(\mathrm{Bi}\), and \(\mathrm{Bi}'\) have infinitely many zeros, all contained in \((-\infty,0)\). In particular, \(\mathrm{Ai}(x)>0\) and \(\mathrm{Bi}(x)>0\) for \(x\ge 0\). In the complex plane, all zeros of \(\mathrm{Ai}\) and \(\mathrm{Ai}'\) lie on the negative real axis, whereas \(\mathrm{Bi}\) and \(\mathrm{Bi}'\) have additional infinite families of zeros in the sectors \[\frac{\pi}{3}<|\arg z|<\frac{\pi}{2},\] symmetric with respect to the real axis.

Every zero of a nontrivial solution of [eq:Airy] is simple: if \(y(z_0)=y'(z_0)=0\), then by uniqueness of the Cauchy problem one must have \(y\equiv 0\). Moreover, \(\mathrm{Ai}\) and \(\mathrm{Bi}\) share no common zeros, since their Wronskian is constant and nonzero (see [wroskiano]). On \((-\infty,0)\) the zeros of linearly independent solutions interlace: between two consecutive zeros of \(\mathrm{Ai}\) there is exactly one zero of \(\mathrm{Bi}\), and between two consecutive zeros of \(\mathrm{Ai}\) there is exactly one zero of \(\mathrm{Ai}'\) (and analogously for \(\mathrm{Bi}\) and \(\mathrm{Bi}'\)).

Asymptotic approximations of the real zeros are classical. As \(k\to\infty\), \[a_k \sim \Big(\tfrac{3\pi}{2}\big(k-\tfrac{1}{4}\big)\Big)^{2/3},\qquad a'_k \sim \Big(\tfrac{3\pi}{2}\big(k-\tfrac{3}{4}\big)\Big)^{2/3},\] \[b_k \sim \Big(\tfrac{3\pi}{2}\big(k-\tfrac{3}{4}\big)\Big)^{2/3},\qquad b'_k \sim \Big(\tfrac{3\pi}{2}\big(k-\tfrac{1}{4}\big)\Big)^{2/3}.\] More precisely, \[a_k = -T\!\Big(\tfrac{3\pi}{8}(4k-1)\Big),\qquad a'_k = -U\!\Big(\tfrac{3\pi}{8}(4k-3)\Big),\] \[b_k = -T\!\Big(\tfrac{3\pi}{8}(4k-3)\Big),\qquad b'_k = -U\!\Big(\tfrac{3\pi}{8}(4k-1)\Big),\] where \(T\) and \(U\) admit complete asymptotic expansions in descending powers and satisfy \(T(x),U(x)\sim x^{2/3}\) as \(x\to\infty\).

A counting function description follows from the oscillatory expansion of \(\mathrm{Ai}(-X)\) as \(X\to\infty\). If \(N_{\mathrm{Ai}}(X)=\#\{k:\ a_k\le X\}\), then \[N_{\mathrm{Ai}}(X)=\frac{1}{\pi}\Big(\frac{2}{3}X^{3/2}+\frac{\pi}{4}\Big)+o(1) =\frac{2}{3\pi}X^{3/2}+\frac{1}{4}+o(1),\] with the shift \(1/4\mapsto 3/4\) for \(N_{\mathrm{Bi}}(X)\).

For reference, the first few real zeros are approximately \[-a_1\approx -2.33811,\quad -a_2\approx -4.08795,\quad -a_3\approx -5.52056,\] \[-b_1\approx -1.17371,\quad -b_2\approx -3.27109,\quad -b_3\approx -4.83074.\] Figure 4 displays \(\mathrm{Ai}(x)\) on the real line together with its zeros at \(x=-a_k\) and \(x=-a'_k\). The compact table embedded at the top of the panel lists the corresponding numerical values used throughout the paper.

Real zeros of \(\operatorname{Ai}(x)\) (filled markers at \(x=-a_k\)) and of \(\operatorname{Ai}'(x)\) (hollow markers at \(x=-a'_k\)). The clean table at the top of the panel reports the numerical values of the first zeros.

Results

Euler–Bernoulli beam supported under its own weight

We consider a homogeneous slender column of length \(L\), oriented vertically, with the lower end clamped at \(x=0\) and the upper end free at \(x=L\). The variables and parameters are defined as follows:

Under the Euler–Bernoulli assumptions—(i) cross sections remain plane and orthogonal to the deformed axis, (ii) deflections and rotations are small (geometric linearization, see  ), (iii) the material follows Hooke’s law, and (iv) arc length \(s\) coincides with the abscissa \(x\) to first order—the axial compressive force at section \(x\) due to self–weight is \[\label{eq:axial-selfweight} N(x)=\int_x^L \delta g\,ds \;\simeq\; \delta g\,(L-x),\] which decreases linearly from the base to the top. The curvature satisfies \(\kappa(x)=d\theta/ds\simeq \theta'(x)\), and the bending moment is \(M(x)=E I\,\kappa(x)\simeq E I\,\theta'(x)\). Linear equilibrium of a differential element, under an axial force \(N(x)\) applied at angle \(\theta\), yields (to first order, using \(\sin\theta\simeq\theta\)) the balance equation \[\label{eq:beam-angle-equil} E I\,\theta''(x)+ N(x)\,\theta(x)=0.\] Substituting [eq:beam-angle-equil] gives the governing model \[\label{eq:columna-ode} E I\,\theta''(x)+ \delta g\,(L-x)\,\theta(x)=0,\qquad 0<x<L,\] with boundary conditions \[\label{eq:bc-phys} \begin{aligned} \theta(0)&=0 \quad &&\text{(clamped end: zero rotation)},\\ M(L)=E I\,\theta'(L)&=0 \quad &&\text{(free end: zero bending moment)}. \end{aligned}\] Introducing the scaling parameter \[\label{eq:a-parameter} a^2:=\frac{\delta g}{E I}\,,\] equation [eq:beam-angle-equil] takes the form \[\theta''(x)+ a^2(L-x)\,\theta(x)=0.\]

Uniform column of length \(L\) under its own weight: sign conventions and variables. The axial force \(N(x)=\delta g(L-x)\) induces flexural instability beyond a critical length, see .

Boundary value problem

The angular deflection \(\theta(x)\) at a point \(P(x)\) satisfies the boundary value problem \[\label{eq:columna-bvp} \begin{aligned} E I\,\theta''(x) + \delta g\,(L-x)\,\theta(x) &= 0,\\ \theta(0) &= 0,\\ \theta'(L) &= 0, \end{aligned}\] where \(E\) is Young’s modulus, \(I\) the moment of inertia of the cross section, \(\delta\) the constant linear density, \(g\) the gravitational acceleration, and \(x\) the distance from the base. Nontrivial solutions exist only for certain values of \(L\), corresponding to the critical buckling lengths.

(a) Change of variable and general solution in terms of Bessel functions

From [eq:columna-bvp], define \[\label{eq:change-vars-21} t := L - x \in (0,L),\qquad \theta(x) =: \phi(t).\] Since \(\tfrac{d}{dx}=-\tfrac{d}{dt}\) and \(\tfrac{d^2}{dx^2}=\tfrac{d^2}{dt^2}\), we obtain \[\label{eq:phi-ode-21} E I\,\phi''(t) + \delta g\,t\,\phi(t) = 0,\qquad \phi(L)=0,\qquad \phi'(0)=0.\] With the scaling parameter [eq:a-parameter], \[\label{eq:phi-airy-21} \phi''(t) + a^2\,t\,\phi(t)=0,\qquad a^2=\frac{\delta g}{E I}.\] Using the canonical mapping introduced in Section 2.1, see Proposition 1, for \(t>0\) we set \[\label{eq:map-21} z=\frac{2}{3}\,a\,t^{3/2},\qquad \phi(t)= t^{1/2} u(z).\] Equation [eq:phi-airy-21] then reduces to the Bessel equation [eq:bessel-clasica] of order \(\nu=\pm\frac{1}{3}\). Reversing the transformation [eq:map-21] yields the consistent basis \[\label{eq:basis-21} \phi(t) = t^{1/2}\!\left[\,C_1\,J_{-1/3}\!\Bigl(\dfrac{2}{3}a\,t^{3/2}\Bigr)\;+\; C_2\,J_{1/3}\!\Bigl(\dfrac{2}{3}a\,t^{3/2}\Bigr)\right].\]

Remark 2 (Equivalence between \(I_{\pm1/3}\)\(J_{\pm1/3}\) via rescaling and analytic continuation). Let \(a>0\). Consider the solution for \(x>0\) \[y(x)=\sqrt{x}\Big\{C_1 I_{1/3}\!\Big(\tfrac{2}{3}x^{3/2}\Big) +C_2 I_{-1/3}\!\Big(\tfrac{2}{3}x^{3/2}\Big)\Big\}.\] Rescaling \(x=a^{2/3}t\) gives \(\xi=\tfrac{2}{3}x^{3/2}=\tfrac{2}{3}a\,t^{3/2}=:z\). If, in addition, we continue into the oscillatory region by taking \(x=-a^{2/3}t\) with \(t>0\), then \(\xi=\tfrac{2}{3}(-a^{2/3}t)^{3/2}=-iz\) and \[I_\nu(-iz)=i^{-\nu}J_\nu(z),\qquad \nu\in\mathbb{C}.\] Since \(\sqrt{x}=\sqrt{-a^{2/3}t}=e^{i\pi/2}a^{1/3}\sqrt{t}\), absorbing phase factors into the constants \(\tilde C_1,\tilde C_2\) yields \[y(-a^{2/3}t)=\sqrt{t}\Big\{\tilde C_1\,J_{-1/3}(z)+\tilde C_2\,J_{1/3}(z)\Big\}, \qquad z=\tfrac{2}{3}a\,t^{3/2},\ \ t>0,\] which coincides with the consistent basis used in [eq:basis-21]. In particular, the representations in terms of \(I_{\pm1/3}\) (exponential region) and \(J_{\pm1/3}\) (oscillatory region) are equivalent under the rescaling \(x=a^{2/3}t\) and analytic continuation across the negative real axis (principal branch choice).

(b) Boundary conditions and critical lengths

As \(t\to 0^+\) (free end), expansions of Bessel functions give \(J_{-1/3}(z)\sim \tfrac{1}{\Gamma(2/3)}(\tfrac{z}{2})^{-1/3}\) and \(J_{1/3}(z)\sim \tfrac{1}{\Gamma(4/3)}(\tfrac{z}{2})^{1/3}\). With \(z=\tfrac{2}{3}a\,t^{3/2}\), [eq:basis-21] becomes \[\phi(t)\sim C_1\,\tfrac{(3/2)^{1/3}}{\Gamma(2/3)a^{1/3}}+ C_2\,\tfrac{(a/6)^{1/3}}{\Gamma(4/3)}\,t+\cdots,\] and the free-end condition \(\phi'(0)=0\) enforces \(C_2=0\). The clamped condition \(\phi(L)=0\) then leads to the spectral equation \[\label{eq:eig-eqn-21} J_{-1/3}\!\Bigl(\dfrac{2}{3}a\,L^{3/2}\Bigr)=0.\] Denoting by \(j_{-1/3,k}\) the positive zeros of \(J_{-1/3}\), the family of critical lengths is \[\label{eq:Lcrit-family-21} L_k=\left(\frac{3\,j_{-1/3,k}}{2\,a}\right)^{2/3},\qquad k=1,2,\dots\] with associated eigenmodes \[\label{eq:theta-mode-21} \theta_k(x)=\phi_k(L-x)= (L_k-x)^{1/2}\,J_{-1/3}\!\Bigl(\dfrac{2}{3}a\,(L_k-x)^{3/2}\Bigr).\] The smallest length \(L_1\) determines the critical threshold \(L_{\mathrm{crit}}\).

(c) Numerical evaluation

For a solid steel rod of radius \(r=0.05\,\mathrm{in}\) with \(A=\pi r^2\), \(I=\tfrac{1}{4}\pi r^4\), \(\delta g=0.28\,A\,\mathrm{lb/in}\), and \(E=2.6\times 10^7\,\mathrm{lb/in}^2\), the parameters are \[A=0.0078539816\;\mathrm{in}^2,\qquad I=0.0000049087\;\mathrm{in}^4.\] From [eq:a-parameter], \[a^2=\frac{\delta g}{E I}\approx 0.0000172308\ \mathrm{in}^{-3},\qquad a\approx 0.0041509962\ \mathrm{in}^{-3/2}.\] Using the first zero \(j_{-1/3,1}\approx 1.8663508589\) of \(J_{-1/3}\), equation [eq:Lcrit-family-21] yields \[\label{eq:Lcrit-num-21} L_{\mathrm{crit}}\approx 76.905\ \mathrm{in}\;\approx\; 6.41\ \mathrm{ft}\;\approx\; 1.953\ \mathrm{m}.\] This value sets the critical instability threshold under self–weight for the idealized column with the given parameters.

Quantum models: linear potential and the quantum bouncer

Several one–dimensional quantum mechanical models reduce, after a suitable rescaling, to Airy-type equations. A particularly relevant case is the linear potential, which provides a first-order approximation of a uniform external field such as gravity near the Earth’s surface. The corresponding spectral problem can thus be formulated in terms of Airy functions, linking directly with the analytical results of Section 2 on power series, Bessel connections, and integral representations (see ).

Consider the stationary Schrödinger equation \[\label{eq:qm-schrodinger-intro} -\frac{\hbar^2}{2M}\,\psi''(z)+V(z)\,\psi(z)=E\,\psi(z),\] with \(V(z)=Fz\) in a region where the field is uniform (\(F>0\) constant). Introducing the characteristic length and the dimensionless variables \[\label{eq:qm-escala-intro} \ell_F:=\Bigl(\frac{\hbar^2}{2M F}\Bigr)^{\!1/3},\qquad \zeta:=\frac{z}{\ell_F}-\varepsilon,\qquad \varepsilon:=\frac{E}{F\,\ell_F},\] and using \(d/dz=(1/\ell_F)\,d/d\zeta\), equation [eq:qm-schrodinger-intro] becomes [eq:Airy] namely, the homogeneous Airy equation already studied in Section 2. Its fundamental solutions are \(\{\mathrm{Ai}(\zeta),\mathrm{Bi}(\zeta)\}\), and the physical selection is dictated by their asymptotic properties: \(\mathrm{Ai}(\zeta)\) decays exponentially for \(\zeta\to+\infty\), while \(\mathrm{Bi}(\zeta)\) diverges.

A paradigmatic realization of this structure is the so-called quantum bouncer, where \(F=Mg\) and the potential includes a rigid wall at \(z=0\), \[\label{eq:qb-pot} V(z)=\begin{cases} Mg\,z,& z>0,\\[2pt] +\infty,& z\le 0. \end{cases}\] The stationary Schrödinger equation is then \[\label{eq:qb-schrodinger} -\frac{\hbar^2}{2M}\,\psi''(z)+Mg\,z\,\psi(z)=E\,\psi(z),\qquad z>0,\] with boundary conditions \[\label{eq:qb-BCs} \psi(0)=0,\qquad \psi(z)\xrightarrow[z\to+\infty]{}0.\] Introducing the length scale \[\label{eq:qb-escala} \ell:=\Bigl(\frac{\hbar^2}{2M^2g}\Bigr)^{\!1/3},\qquad \zeta:=\frac{z}{\ell}-\varepsilon,\qquad \varepsilon:=\frac{E}{Mg\,\ell},\] and rescaling derivatives accordingly, equation [eq:qb-schrodinger] reduces to [eq:Airy], again the homogeneous Airy equation. Square-integrability forces the rejection of \(\mathrm{Bi}\), leaving \[\psi(\zeta)=C\,\mathrm{Ai}(\zeta),\qquad C\in\mathbb{C}.\] The rigid wall at \(z=0\) corresponds to \(\zeta(0)=-\varepsilon\), and the boundary condition [eq:qb-BCs] requires \[\label{eq:qb-quant-cond} \mathrm{Ai}(-\varepsilon)=0.\] Let \(\{-a_n\}_{n\ge 1}\) denote the negative zeros of \(\mathrm{Ai}\), ordered as \(0<a_1<a_2<\cdots\), see Section 2.3. Then \(\varepsilon=a_n\) and the energy spectrum is \[\label{eq:qb-energies} E_n=Mg\,\ell\,a_n=\Bigl(\tfrac{1}{2}M g^2\hbar^2\Bigr)^{\!1/3}\,a_n=:E_0\,a_n,\qquad n=1,2,\dots,\] with associated (unnormalized) eigenfunctions \[\label{eq:qb-modes} \psi_n(z)=\mathrm{Ai}\!\Bigl(\tfrac{z}{\ell}-a_n\Bigr),\qquad z\ge 0.\]

The integral representation of \(\mathrm{Ai}\) established in Section 2.2 justifies its exponential decay for \(\zeta\to+\infty\). The quantization condition [eq:qb-quant-cond] can be interpreted as a “phase selection’’ of the oscillatory integral (steepest–descent contours): only for \(\varepsilon=a_n\) does \(\mathrm{Ai}(-\varepsilon)\) vanish, thus enforcing the boundary condition at \(z=0\). Moreover, the characteristic length \(\ell\) sets the decay scale near the wall, while \(E_0=(\tfrac{1}{2}M g^2\hbar^2)^{1/3}\) fixes the spacing between consecutive energy levels.

Finally, normalization is achieved by expressing the constants \(C_n\) in terms of \(\mathrm{Ai}'(-a_n)\) via the identity \[\int_{-a_n}^{+\infty}\mathrm{Ai}^2(\zeta)\,d\zeta=\mathrm{Ai}'^2(-a_n).\] With the change \(z=\ell(\zeta+a_n)\) one obtains \[C_n=\bigl(\ell\,\mathrm{Ai}'^2(-a_n)\bigr)^{-1/2},\] so that the functions [eq:qb-modes] form an orthonormal system in \(L^2(0,\infty)\). In practice, the zeros \(a_n\) and derivatives \(\mathrm{Ai}'(-a_n)\) are computed numerically and inserted in [eq:qb-modes].

Particle in a uniform electric field

We consider a particle of mass \(m\) and charge \(q\) moving in a region where the electric field is uniform and directed along the \(+\hat{x}\) axis. The corresponding scalar potential is taken as \[\label{eq:e-uniform-pot} V(x)=q\,E\,x,\qquad E>0,\] so that the potential energy increases linearly with \(x\). The stationary Schrödinger equation reads \[\label{eq:e-uniform-sch} -\frac{\hbar^2}{2m}\,\psi''(x)+qE\,x\,\psi(x)=\mathcal{E}\,\psi(x),\] where \(\mathcal{E}\) denotes the energy eigenvalue. If \(qE<0\), the substitution \(x\mapsto -x\) restores the case [eq:e-uniform-pot] without loss of generality.

Introducing the characteristic length scale and the dimensionless variable \[\label{eq:e-uniform-scale} \ell_E:=\Bigl(\frac{\hbar^2}{2m\,qE}\Bigr)^{\!1/3},\qquad \zeta:=\frac{x}{\ell_E}-\varepsilon,\qquad \varepsilon:=\frac{\mathcal{E}}{qE\,\ell_E},\] we have \(\tfrac{d}{dx}=(1/\ell_E)\tfrac{d}{d\zeta}\) and \(\tfrac{d^2}{dx^2}=(1/\ell_E^2)\tfrac{d^2}{d\zeta^2}\), and equation [eq:e-uniform-sch] reduces to [eq:Airy] , consistent with Section 2.

The general solution of [eq:e-uniform-sch] is a linear combination of Airy functions, \[\label{eq:e-uniform-sol} \psi(x)=A\,\mathrm{Ai}\!\Bigl(\dfrac{x}{\ell_E}-\varepsilon\Bigr)+B\,\mathrm{Bi}\!\Bigl(\dfrac{x}{\ell_E}-\varepsilon\Bigr),\qquad A,B\in\mathbb{C}.\] For \(x\to+\infty\) (equivalently \(\zeta\to+\infty\)) the function \(\mathrm{Bi}\) grows exponentially, while \(\mathrm{Ai}\) decays exponentially, as follows from the integral representation and asymptotics of Section 2.2. Consequently, for square-integrable states or physically acceptable radiative boundary conditions one imposes \(B=0\).

When the system is confined to \(x\geq 0\) by a rigid wall at \(x=0\), corresponding to an electrostatic triangular well, the boundary condition \[\label{eq:e-uniform-wall} \psi(0)=0\] leads to the spectral equation \[\label{eq:e-uniform-quant} \mathrm{Ai}(-\varepsilon)=0.\] Denoting the negative zeros of \(\mathrm{Ai}\) by \(\{-a_n\}_{n\ge 1}\) with \(0<a_1<a_2<\cdots\), the quantization condition \(\varepsilon=a_n\) yields the discrete spectrum \[\label{eq:e-uniform-energies} \mathcal{E}_n=qE\,\ell_E\,a_n=\Bigl(\tfrac{1}{2}m\,(qE)^2\hbar^2\Bigr)^{\!1/3}a_n=:\mathcal{E}_0\,a_n,\qquad n=1,2,\dots,\] together with the associated eigenfunctions \[\label{eq:e-uniform-modes} \psi_n(x)=\mathrm{Ai}\!\Bigl(\dfrac{x}{\ell_E}-a_n\Bigr),\qquad x\ge 0,\] where the normalization constant is determined by \(\mathrm{Ai}'(-a_n)\), in complete analogy with the gravitational “quantum bouncer” of Section 2.3 (cf. [eq:qb-modes]).

From a physical viewpoint, the model [eq:e-uniform-sch] is the electrostatic analogue of the quantum bouncer with the identification \(qE\equiv Mg\). In semiconductor physics, the quantization condition [eq:e-uniform-quant] describes bound states in triangular quantum wells (depletion regions or heterostructure interfaces), while the general solution [eq:e-uniform-sol] governs tunneling through linearized barriers, as in Fowler–Nordheim emission. Moreover, the integral representation of Section 2.2 provides a rigorous foundation for the asymptotic selection of [eq:e-uniform-sol] and for the construction of Green’s functions in forced problems with linear potentials (see ).

The Schrödinger equation for a charged particle in a uniform electric field leads to solutions expressed in terms of Airy functions. This mathematical structure is not only of theoretical interest but has also been experimentally realized in the context of electron Airy beams (see ).

Discussion

The results obtained here should be interpreted at the level of analytical structure rather than as a collection of isolated formulas. Although the mathematical ingredients developed in Section 2 are classical, their role in this paper is not merely expository. Instead, they are organized into a coherent derivational framework through which a common Airy-function mechanism becomes explicit across several canonical models with linear-potential character.

From this standpoint, the contribution of the manuscript is not to reintroduce Airy functions as objects of classical special-function theory, but to demonstrate how their series construction, Bessel reductions, contour-integral representations, and zero structure act in concert in the derivation of physically meaningful spectral and critical conditions. At this level, the analysis acquires structural significance: the same mathematical framework governs admissible branches, asymptotic selection, and characteristic scales in continuum and quantum settings that are often treated independently in the literature.

This perspective also clarifies the position of the paper with respect to earlier references. Standard sources such as , , and provide the classical background in special-function theory, while , , and emphasize the role of Airy functions in linear-potential quantum mechanics and turning-point analysis. Semiconductor applications are addressed in and , whereas the quantum bouncer appears as a specific spectral model in the physics literature. What the present work contributes is not a bibliographic synthesis of these directions, but an explicit analytical scheme that places them within a unified framework and shows how the relevant spectral information emerges from a common Airy-function structure.

In particular, the paper establishes a direct correspondence between mathematical representation and physical interpretation. The power-series construction determines the canonical solutions and their normalizations; the Bessel reductions separate the exponential and oscillatory regimes; the contour formulation provides analytic continuation and asymptotic control; and the zero distribution of \(\mathrm{Ai}\) governs quantization and criticality in the models under consideration. This integrated viewpoint constitutes the main analytical contribution of the manuscript.

Accordingly, the paper should be read neither as a broad review nor as a purely pedagogical account. Its aim is to provide a rigorous and self-contained analytical treatment of a selected class of models for which Airy functions are structurally decisive. In this sense, its value lies in making explicit, within a single derivational setting, the common mechanism by which special-function theory governs representative problems in continuum mechanics and one-dimensional quantum mechanics.

Conclusions

This paper establishes a coherent analytical framework for Airy functions and demonstrates how it governs a selected class of canonical models in continuum and quantum mechanics. Starting from the Airy equation \(y''-x\,y=0\) [eq:Airy], we derived the power-series construction, the corresponding normalizations at the origin, the reductions to ordinary and modified Bessel equations, and the contour-integral representations associated with the canonical solutions \(\mathrm{Ai}\) and \(\mathrm{Bi}\). These complementary formulations provide a unified description of the oscillatory and exponential regimes, together with the asymptotic structure required for spectral and boundary-value analysis.

The main outcome is that this Airy framework is not merely formal, but structurally decisive in the models considered. For the Euler–Bernoulli column under self-weight, the reduction to Bessel form yields the spectral condition [eq:eig-eqn-21], from which the critical buckling lengths [eq:Lcrit-family-21] and associated modes [eq:theta-mode-21] follow. For the quantum bouncer and the particle in a uniform electric field, the corresponding rescalings reduce Schrödinger’s equation to Airy form, and the admissible spectrum is determined by the zeros of \(\mathrm{Ai}\), leading to the quantization conditions [eq:qb-quant-cond] and [eq:e-uniform-quant], the spectra [eq:qb-energies] and [eq:e-uniform-energies], and the eigenmodes [eq:qb-modes] and [eq:e-uniform-modes]. In all cases, the relevant physical scales emerge directly from the Airy reduction.

Taken together, these results show that Airy functions provide an exact analytical bridge between special-function theory and representative models with linear-potential structure. Their series, Bessel, and contour representations are mutually consistent and analytically effective, allowing one to derive spectra, critical thresholds, and mode shapes within a single framework. This also suggests a natural extension toward more general variable-coefficient, perturbed, or mixed-boundary problems, where the Airy structure is expected to remain a fundamental organizing principle.

Declarations

Authors’ contributions

All authors contributed to the preparation of this manuscript. Juan Toribio Milane led the development of the quantum and applied physical models. José Angel Gómez Hernández derived the Bessel connections and the contour-integral representations. Francis Leandro Álvarez Paulino and Manuel Leonardo Reyes Cordero developed and analyzed the solutions of the Airy equation. All authors contributed to the discussion of the results and approved the final version of the manuscript.

How to cite this article

Milane Toribio, J., Gómez Hernández, J. A., Álvarez Paulino, F. L., & Reyes Cordero, M. L. (2025). Models of Beams and Quantum Systems with Linear Potentials via Airy Functions. Revista Alma Mater, [volume(issue)], [pp–pp]. https://doi.org/[DOI]

Funding

The research of Juan Toribio Milane was partially supported by the Fondo Nacional de Innovación y Desarrollo Científico y Tecnológico (FONDOCYT), Dominican Republic, under grant number 2024-2-1D2-0791.

Competing interests

The authors declare that they have no competing interests.

Generative AI use

No generative AI tools were used in the preparation of this manuscript.

99

Cacciapuoti, Claudio, et al. (2025). The Born-Oppenheimer approximation for a 1D 2+1 particle system with zero-range interactions. Journal of Mathematical Physics. https://arxiv.org/abs/2506.21457

Wine, N., Achtymichuk, J., & Marsiglio, F. (2025). The modified Airy function approximation applied to the double-well potential. AIP Advances, 15(3), 035330. https://doi.org/10.1063/5.0241523

Eberhardt, A., & Hui, L. (2025). Caustic fringes for wave dark matter. https://arxiv.org/pdf/2506.02400

Davies, J. H. (1998). The Physics of Low-Dimensional Semiconductors: An Introduction. Cambridge University Press. (see chapter: Airy Functions: Triangular Well)

Landau, L. D., & Lifshitz, E. M. (1977). Quantum Mechanics: Non-Relativistic Theory (3rd ed., revised and enlarged). Pergamon Press.

Zakrzewski, M., & Żołądek, H. (2025). Variations on hypergeometric functions. https://doi.org/10.48550/arXiv.2501.08310

Abramochkin, E. G., & Razueva, E. V. (2018). Higher Derivatives of Airy Functions and of their Products. SIGMA, 14, 042. https://doi.org/10.3842/SIGMA.2018.042

Siviloglou, G. A., Broky, J., Dogariu, A., & Christodoulides, D. N. (2007). Observation of accelerating Airy beams. Physical Review Letters, 99, 213901. https://doi.org/10.1103/PhysRevLett.99.213901

VolochBloch, N., Lereah, Y., Lilach, Y., Gover, A., & Arie, A. (2013). Generation of electron Airy beams. Nature, 494(7437), 331–335. https://doi.org/10.1038/nature11840

Olver, F. W. J., Lozier, D. W., Boisvert, R. F., & Clark, C. W. (2010). NIST Handbook of Mathematical Functions. Cambridge University Press. ISBN 978-0521192255.

Vallée, O., & Soares, M. (2010). Airy Functions and Applications to Physics (2nd ed.). Imperial College Press / World Scientific. ISBN 978-1848165489.

Gradshteyn, I. S., & Ryzhik, I. M. (2014). Table of Integrals, Series, and Products (8th ed.; D. Zwillinger & V. H. Moll, Eds.). Academic Press. ISBN 978-0-12-384933-5.

Griffiths, D. J., & Schroeter, D. F. (2018). Introduction to Quantum Mechanics (3rd ed.). Cambridge University Press. ISBN 978-1107189638.

Murphy, E. L., & Good, R. H. Jr. (1964). WKB Connection Formulas. Journal of Mathematics and Physics (now Studies in Applied Mathematics), 43(1–4), 251–254. https://doi.org/10.1002/sapm1964431251

Villavicencio, F. L. (2018). Estudio de un modelo teórico para el doblamiento de bandas en un semiconductor [Study of Band Bending in a Semiconductor]. Bachelor’s thesis (Licenciatura en Física), https://www.facet.unt.edu.ar/licfisica/wp-content/uploads/sites/57/2018/12/Tesis-FacundoVillavicencio-sept2018.pdf (accessed March 2025).

Agarwal, R. P., & O’Regan, D. (2009). Ordinary and Partial Differential Equations. Springer.

Abramowitz, M., & Stegun, I. A. (1972). Handbook of Mathematical Functions. Dover.

Ahlfors, L. V. (1979). Complex Analysis (3rd ed.). McGraw–Hill.

Dingle, R. B. (1973). Asymptotic Expansions: Their Derivation and Interpretation. Academic Press.

NIST Digital Library of Mathematical Functions. F. W. J. Olver et al. (eds.). https://dlmf.nist.gov, Release 1.2.4 of 2025-03-15.

GeaBanacloche, J. (1999). A quantum bouncing ball. American Journal of Physics, 67(9), 776–782. https://doi.org/10.1119/1.19124

Rudin, W. (1987). Real and Complex Analysis (3rd ed.). McGraw–Hill. ISBN 0-07-054234-1.

Vankov, A. A. (2009). Quantum bouncer: theory and experiment. https://doi.org/10.48550/arXiv.0906.5138.

Lebedev, N. N. (1972). Special Functions and Their Applications (R. A. Silverman, Trans.). Dover Publications, New York.

Riley, K. F., Hobson, M. P., & Bence, S. J. (2006). Mathematical Methods for Physics and Engineering: A Comprehensive Guide (3rd ed.). Cambridge University Press. ISBN 0521679710

Zill, D. G. (2014). Differential Equations with Modeling Applications (Vol. 1, 10th ed.). Cengage Learning. [In Spanish: Ecuaciones diferenciales con aplicaciones de modelado].