Neumann Series of Bessel Functions (NSBF)

• NSBF are series expansions representing differential equation solutions using Bessel or spherical Bessel functions with coefficients encoding operator or potential information.
• They provide uniform convergence with respect to the spectral parameter, ensuring stable and accurate numerical results in applications like Schrödinger, Sturm–Liouville, and Dirac equations.
• The method leverages transmutation operators, Legendre expansions, and recursive integration schemes to construct reliable approximations and efficient inverse problem solutions.

Neumann series of Bessel functions (NSBF) are representations in which a function, or more specifically a solution of a differential equation, is expanded in Bessel or spherical Bessel functions with coefficients that encode the underlying operator or potential. In the classical sense, an NSBF has the form $$\mathfrak N_\nu(z)\sim \sum_{n=0}^\infty a_n J_{\nu+n}(z)$$ (Micheli, 2017). In the modern transmutation-based literature, the term usually refers to representations of Schrödinger, Sturm–Liouville, perturbed Bessel, and Dirac solutions as series in $j_n(\omega x)$ or related Bessel models whose truncation error is uniform with respect to the spectral parameter on the real axis, and in horizontal strips of the complex plane (Kravchenko et al., 2015, Kravchenko et al., 2016, Kravchenko et al., 2020, Roque et al., 1 Feb 2025). This uniform-in-spectral-parameter behavior is the defining practical feature in the computational branch of the subject.

1. Definition and scope

The materials surveyed here show two recurrent usages of the term. The first is the classical special-function usage, where one studies series directly in Bessel functions of shifted order,

$$\mathfrak N_\nu(z)\sim \sum_{n=0}^\infty a_n J_{\nu+n}(z),$$

and seeks integral representations, convergence criteria, or closed-form sums (Micheli, 2017). In that setting, the coefficient sequence is encoded by the analytic generating function $A(z)=\sum_{n=0}^\infty a_n z^n$, and the NSBF can be written as a single integral against a universal Bessel kernel $K_z^{(\nu)}(\theta)$ (Micheli, 2017).

A second classical branch studies “Neumann-type” series built from products of Bessel functions. The paper on Turán determinants introduces

$$\mathfrak G_{\mu,\nu}^{a,b}(x):=\sum_{n\ge 1}\theta_n J_{\mu+an}(x)J_{\nu+bn}(x),$$

calling these “second kind Neumann type series of Bessel functions of the first kind” (Baricz et al., 2011). In that literature, von Lommel’s formula and related Turán identities generate positive series in squares $J_{\nu+1+2n}^2(x)$, and the main structural issue is the conversion of such sums into double-integral representations (Baricz et al., 2011).

A third, later usage is computational and operator-theoretic. Here NSBF denotes an exact representation of a solution of a differential equation as a free oscillatory term plus a series in spherical Bessel functions with coefficient functions depending only on the spatial variable and the coefficients of the equation (Kravchenko et al., 2015, Kravchenko et al., 2016, Kravchenko et al., 2016, Kravchenko et al., 2020, Roque et al., 1 Feb 2025). This is the usage that dominates current spectral and inverse-problem applications.

The term is not synonymous with every infinite expansion involving Bessel functions. One paper included here explicitly states that its formulas for $J_n$, derived as Fourier series of auxiliary functions and expressed through elementary trigonometric functions of $\sqrt{x^2+(k\pi)^2}$, are not Neumann series of Bessel functions in the classical sense (Andrusyk, 2012). That distinction is central: NSBF are expansions in a Bessel basis, not merely expansions of Bessel functions by some other oscillatory basis.

2. Analytic framework: transmutations, formal powers, and Legendre coefficients

The modern NSBF method is built on transmutation operators. For the one-dimensional Schrödinger equation

$-y''+q(x)y=\omega^2 y,$

the transmutation operator is

$j_n(\omega x)$0

with kernel $j_n(\omega x)$1 independent of $j_n(\omega x)$2, and it maps free solutions of $j_n(\omega x)$3 to solutions of the perturbed equation (Kravchenko et al., 2015). This $j_n(\omega x)$4-independence is the structural reason that truncated NSBF approximations are uniform in the spectral parameter.

The kernel is accessed through formal powers rather than by solving for $j_n(\omega x)$5 directly. A nonvanishing solution $j_n(\omega x)$6 of $j_n(\omega x)$7, normalized by $j_n(\omega x)$8, generates recursive integrals $j_n(\omega x)$9, $$\mathfrak N_\nu(z)\sim \sum_{n=0}^\infty a_n J_{\nu+n}(z),$$0, and the formal powers $$\mathfrak N_\nu(z)\sim \sum_{n=0}^\infty a_n J_{\nu+n}(z),$$1, $$\mathfrak N_\nu(z)\sim \sum_{n=0}^\infty a_n J_{\nu+n}(z),$$2 used throughout SPPS and NSBF theory (Kravchenko et al., 2015). The mapping property

$$\mathfrak N_\nu(z)\sim \sum_{n=0}^\infty a_n J_{\nu+n}(z),$$3

links these recursive integrals to the transmutation operator (Kravchenko et al., 2015).

The decisive analytic step is the Fourier–Legendre expansion of the transmutation kernel,

$$\mathfrak N_\nu(z)\sim \sum_{n=0}^\infty a_n J_{\nu+n}(z),$$4

with coefficients

$$\mathfrak N_\nu(z)\sim \sum_{n=0}^\infty a_n J_{\nu+n}(z),$$5

and an analogous expansion for $$\mathfrak N_\nu(z)\sim \sum_{n=0}^\infty a_n J_{\nu+n}(z),$$6 with coefficients $$\mathfrak N_\nu(z)\sim \sum_{n=0}^\infty a_n J_{\nu+n}(z),$$7 (Kravchenko et al., 2015). This Legendre expansion is what converts the transmutation formula into an NSBF, because integrals of $$\mathfrak N_\nu(z)\sim \sum_{n=0}^\infty a_n J_{\nu+n}(z),$$8 against $$\mathfrak N_\nu(z)\sim \sum_{n=0}^\infty a_n J_{\nu+n}(z),$$9, $A(z)=\sum_{n=0}^\infty a_n z^n$0, or $A(z)=\sum_{n=0}^\infty a_n z^n$1 are expressible through spherical Bessel functions (Kravchenko et al., 2015, Kravchenko et al., 2016).

The same pattern survives under Liouville transformation for general Sturm–Liouville equations. Writing

$A(z)=\sum_{n=0}^\infty a_n z^n$2

one transforms the Sturm–Liouville equation to a Schrödinger equation, constructs NSBF there, and then pulls the representation back to the original variables (Kravchenko et al., 2016). In the Dirac case, the transmutation operator becomes matrix-valued; the one-dimensional Dirac paper expands its kernel as

$A(z)=\sum_{n=0}^\infty a_n z^n$3

and obtains matrix NSBF coefficients $A(z)=\sum_{n=0}^\infty a_n z^n$4 from this expansion (Roque et al., 1 Feb 2025).

3. Canonical representations

For the one-dimensional Schrödinger equation, the fundamental solutions normalized at the left endpoint are represented as

$A(z)=\sum_{n=0}^\infty a_n z^n$5

$A(z)=\sum_{n=0}^\infty a_n z^n$6

with corresponding derivative series involving $A(z)=\sum_{n=0}^\infty a_n z^n$7 (Kravchenko et al., 2015). This is the prototype from which much of the subsequent literature proceeds.

A refined large-$A(z)=\sum_{n=0}^\infty a_n z^n$8 variant isolates the first asymptotic terms explicitly and represents only the remainder by an NSBF: $A(z)=\sum_{n=0}^\infty a_n z^n$9 with $K_z^{(\nu)}(\theta)$0, $K_z^{(\nu)}(\theta)$1, and $K_z^{(\nu)}(\theta)$2 (Kravchenko et al., 2017). The leading asymptotic structure is thus explicit, while the remainder retains the constructive NSBF form.

For singular equations of perturbed Bessel type,

$K_z^{(\nu)}(\theta)$3

the regular solution has the representation

$K_z^{(\nu)}(\theta)$4

and its derivative has an analogous series with coefficients $K_z^{(\nu)}(\theta)$5 (Kravchenko et al., 2016). Here the unperturbed Bessel term is exact when $K_z^{(\nu)}(\theta)$6, and the NSBF correction records the perturbation.

For general Sturm–Liouville equations, the Liouville pullback yields

$K_z^{(\nu)}(\theta)$7

$K_z^{(\nu)}(\theta)$8

together with derivative expansions involving coefficients $K_z^{(\nu)}(\theta)$9 (Kravchenko et al., 2016). The coefficient functions are Liouville pullbacks of the Schrödinger-side $$\mathfrak G_{\mu,\nu}^{a,b}(x):=\sum_{n\ge 1}\theta_n J_{\mu+an}(x)J_{\nu+bn}(x),$$0 and $$\mathfrak G_{\mu,\nu}^{a,b}(x):=\sum_{n\ge 1}\theta_n J_{\mu+an}(x)J_{\nu+bn}(x),$$1.

The radial Dirac system of special form is reduced to a pair of second-order perturbed Bessel equations by decoupling, with

$$\mathfrak G_{\mu,\nu}^{a,b}(x):=\sum_{n\ge 1}\theta_n J_{\mu+an}(x)J_{\nu+bn}(x),$$2

as the scalar spectral quantity (Kravchenko et al., 2020). The regular solution then has the NSBF form

$$\mathfrak G_{\mu,\nu}^{a,b}(x):=\sum_{n\ge 1}\theta_n J_{\mu+an}(x)J_{\nu+bn}(x),$$3

$$\mathfrak G_{\mu,\nu}^{a,b}(x):=\sum_{n\ge 1}\theta_n J_{\mu+an}(x)J_{\nu+bn}(x),$$4

with explicit zero-energy seed solutions $$\mathfrak G_{\mu,\nu}^{a,b}(x):=\sum_{n\ge 1}\theta_n J_{\mu+an}(x)J_{\nu+bn}(x),$$5 and $$\mathfrak G_{\mu,\nu}^{a,b}(x):=\sum_{n\ge 1}\theta_n J_{\mu+an}(x)J_{\nu+bn}(x),$$6 and recursively computable coefficients (Kravchenko et al., 2020).

For the one-dimensional Dirac equation, the normalized fundamental matrix satisfies

$$\mathfrak G_{\mu,\nu}^{a,b}(x):=\sum_{n\ge 1}\theta_n J_{\mu+an}(x)J_{\nu+bn}(x),$$7

where $$\mathfrak G_{\mu,\nu}^{a,b}(x):=\sum_{n\ge 1}\theta_n J_{\mu+an}(x)J_{\nu+bn}(x),$$8 is the free fundamental matrix and $$\mathfrak G_{\mu,\nu}^{a,b}(x):=\sum_{n\ge 1}\theta_n J_{\mu+an}(x)J_{\nu+bn}(x),$$9 are Legendre coefficients of the Dirac transmutation kernel (Roque et al., 1 Feb 2025). The NSBF therefore extends from scalar second-order equations to first-order matrix systems without passing through a scalar reduction.

4. Convergence, truncation, and spectral-parameter uniformity

The central analytic property of NSBF in the computational literature is uniform convergence with respect to the spectral parameter. For the Schrödinger equation, the series for $J_{\nu+1+2n}^2(x)$0 and $J_{\nu+1+2n}^2(x)$1 converge uniformly in $J_{\nu+1+2n}^2(x)$2 and uniformly in $J_{\nu+1+2n}^2(x)$3 on compact subsets of $J_{\nu+1+2n}^2(x)$4; the truncated sums satisfy real-axis bounds

$J_{\nu+1+2n}^2(x)$5

and strip bounds depending only on $J_{\nu+1+2n}^2(x)$6, not on $J_{\nu+1+2n}^2(x)$7 (Kravchenko et al., 2015). This is the precise sense in which the approximation is uniform with respect to the spectral parameter.

The asymptotic NSBF for the normalized Schrödinger solution improves this further. Its truncation error satisfies

$J_{\nu+1+2n}^2(x)$8

for real $J_{\nu+1+2n}^2(x)$9, and a corresponding strip estimate in the complex plane (Kravchenko et al., 2017). The high-frequency regime is therefore not merely nondeteriorating; the estimate improves like $J_n$0.

For perturbed Bessel equations, the regular solution obeys

$J_n$1

for all real $J_n$2, and the same strip-uniform pattern holds in $J_n$3 (Kravchenko et al., 2016). That paper also proves coefficient-decay estimates and identifies a saturation phenomenon for noninteger $J_n$4: for noninteger $J_n$5, no matter how smooth $J_n$6 is, the decay exponent cannot exceed $J_n$7 (Kravchenko et al., 2016).

The Sturm–Liouville pullback preserves these uniform-in-$J_n$8 bounds. The truncations $J_n$9, $\sqrt{x^2+(k\pi)^2}$0 satisfy real-axis and strip estimates of the same form, up to the Liouville factor $\sqrt{x^2+(k\pi)^2}$1 (Kravchenko et al., 2016).

The radial Dirac paper makes the same property explicit in first-order form. Its truncated approximants satisfy

$\sqrt{x^2+(k\pi)^2}$2

for all real $\sqrt{x^2+(k\pi)^2}$3, with $\sqrt{x^2+(k\pi)^2}$4 independent of $\sqrt{x^2+(k\pi)^2}$5, and an analogous statement in horizontal strips (Kravchenko et al., 2020). Under additional smoothness assumptions on $\sqrt{x^2+(k\pi)^2}$6, it also proves algebraic decay

$\sqrt{x^2+(k\pi)^2}$7

The one-dimensional Dirac paper states the same hallmark property for its matrix-valued NSBF: the representation is uniformly convergent with respect to the spectral parameter, and this is the reason the method can compute large sets of eigendata with non-deteriorating accuracy (Roque et al., 1 Feb 2025).

5. Computation and numerical evidence

The algorithmic pattern is strikingly stable across papers. One first computes a nonvanishing zero-energy or homogeneous seed solution, such as $\sqrt{x^2+(k\pi)^2}$8 for Schrödinger, $\sqrt{x^2+(k\pi)^2}$9 for perturbed Bessel, $-y''+q(x)y=\omega^2 y,$0 for Sturm–Liouville, $-y''+q(x)y=\omega^2 y,$1 for the radial Dirac system, or $-y''+q(x)y=\omega^2 y,$2 for the one-dimensional Dirac system (Kravchenko et al., 2015, Kravchenko et al., 2016, Kravchenko et al., 2020, Roque et al., 1 Feb 2025). One then generates NSBF coefficients recursively by quadratures, evaluates truncated Bessel sums, and constructs characteristic functions or residuals from these truncated representations.

A major practical theme is the replacement of unstable direct coefficient formulas by recurrent integration schemes. For the Schrödinger equation, the direct $-y''+q(x)y=\omega^2 y,$3 formulas are exact but the paper derives a more stable $-y''+q(x)y=\omega^2 y,$4-recurrence and reports that it “clearly outperformed the direct $-y''+q(x)y=\omega^2 y,$5 formulas” numerically (Kravchenko et al., 2015). For perturbed Bessel equations, the authors state that the direct formulas are not numerically stable and recommend recursive quadratures that use only integrations and previously computed quantities (Kravchenko et al., 2016). The radial Dirac paper emphasizes the same point: the recursions avoid repeated differentiation of numerically computed coefficients, while the only derivatives required are those of the explicit seed solutions $-y''+q(x)y=\omega^2 y,$6 (Kravchenko et al., 2020). In the one-dimensional Dirac case, the coefficient matrices are computed from a Volterra solution operator $-y''+q(x)y=\omega^2 y,$7, again yielding an integration-based recurrence (Roque et al., 1 Feb 2025).

The numerical record reported in these papers is unusually strong. For the one-dimensional Schrödinger equation, the NSBF method computed 500 eigenvalues for the Paine problem with maximum absolute error $-y''+q(x)y=\omega^2 y,$8 in high precision, 500 eigenvalues for $-y''+q(x)y=\omega^2 y,$9 with maximum relative error $j_n(\omega x)$00, and 40 complex eigenvalues for $j_n(\omega x)$01 with relative errors below $j_n(\omega x)$02 (Kravchenko et al., 2015). For perturbed Bessel equations with $j_n(\omega x)$03 and $j_n(\omega x)$04, the 100th eigenvalue was computed with absolute error about $j_n(\omega x)$05, while SPPS gave error about $j_n(\omega x)$06 in the authors’ comparison (Kravchenko et al., 2016). For the Sturm–Liouville example $j_n(\omega x)$07, the first 100 eigenvalues were computed with maximal absolute error $j_n(\omega x)$08 and maximal relative error $j_n(\omega x)$09 (Kravchenko et al., 2016). For the radial Dirac oscillator, the observed absolute error did not grow when $j_n(\omega x)$10 became large, except near the right endpoint where machine precision affected precomputed coefficients, directly confirming the spectral-parameter-independent error estimates (Kravchenko et al., 2020).

The newest transmission-eigenvalue work extends the same computational philosophy to inverse problems. After Liouville transformation of the spherically symmetric transmission eigenvalue problem, the characteristic function is expanded in NSBF, eigenvalues are computed by root finding on a truncated partial sum, the transformed interval length $j_n(\omega x)$11 is recovered from spectral data by an NSBF-based consistency criterion, and the refractive index $j_n(\omega x)$12 is reconstructed from the first NSBF coefficients (Kravchenko et al., 22 Jul 2025). That paper explicitly states that the inverse algorithm requires no a priori assumptions on $j_n(\omega x)$13 or on the sign of the contrast $j_n(\omega x)$14, and also introduces spectrum completion by recovering additional zeros from the truncated NSBF characteristic function (Kravchenko et al., 22 Jul 2025).

6. Extensions, generalizations, and conceptual boundaries

NSBF methods have migrated well beyond scalar boundary-value problems. In quantitative finance, double barrier knock-out option pricing under one-dimensional time-homogeneous diffusions is reduced to a regular Sturm–Liouville problem, and each eigenfunction is represented as

$j_n(\omega x)$15

which then feeds directly into the pricing series $j_n(\omega x)$16 (Kravchenko et al., 2017). In electromagnetics, time-dependent one-dimensional Maxwell propagation through an inhomogeneous layer is reduced to a Vekua equation, the associated transmutation kernels are expanded in Legendre polynomials, and for trigonometric Fourier input signals the resulting field is written as an NSBF in spherical Bessel functions with exact medium-dependent coefficients $j_n(\omega x)$17, $j_n(\omega x)$18 (Khmelnytskaya et al., 2019).

The subject also admits exact summation and symmetry-theoretic variants. The paper on the Laplacian on polygons evaluates arithmetic-progression Bessel series such as

$j_n(\omega x)$19

as finite sums of trigonometric functions, interpreting them as generalized Neumann expansions of polygonal eigenstates (Molinari, 2020). This is not a transmutation-based NSBF method, but it lies in the same classical family of Bessel-basis summation problems.

At the level of abstraction, the classical Bessel–Neumann expansion can itself be generalized. One paper replaces the standard Bessel basis $j_n(\omega x)$20 by a basis $j_n(\omega x)$21 generated by an infinite-dimensional linear ODE with an upper Hessenberg matrix $j_n(\omega x)$22, with the Bessel case recovered from the tridiagonal matrix encoding

$j_n(\omega x)$23

In that framework, the coefficients are computed from

$j_n(\omega x)$24

and the basis functions from $j_n(\omega x)$25 (Koskela et al., 2017). This does not replace NSBF, but it makes precise how the Bessel–Neumann expansion sits inside a broader Hessenberg-generated class.

A persistent boundary line remains important. Some works derive remarkable infinite series involving Bessel functions without producing NSBF in the strict sense. The 2012 paper on integer-order $j_n(\omega x)$26 explicitly states that its formulas are not classical Neumann series of Bessel functions, because the basis functions are elementary combinations of $j_n(\omega x)$27, $j_n(\omega x)$28, and rational functions of $j_n(\omega x)$29, not $j_n(\omega x)$30 or $j_n(\omega x)$31 themselves (Andrusyk, 2012). The distinction clarifies the field: NSBF are not defined by the mere presence of Bessel functions, but by expansion in a Bessel basis with structurally significant coefficient systems.

Taken together, these works portray NSBF as a hybrid subject linking classical Bessel-series analysis, transmutation operators, SPPS formal powers, Fourier–Legendre kernel expansions, and spectral computation. Its most stable identifying feature is not a single canonical formula, but a recurring architecture: a free Bessel or oscillatory model term, coefficient functions determined by the underlying operator, and truncation errors controlled uniformly in the spectral parameter.