Representation of solutions of the one-dimensional Dirac equation in terms of Neumann series of Bessel functions

This paper presents a uniformly convergent representation of solutions to the one-dimensional Dirac equation as Neumann series of Bessel functions with explicitly calculated coefficients, enabling an efficient numerical method for solving initial-value and spectral problems with high accuracy across large sets of eigendata.

The Gist

The Big Picture: Tuning a Radio in a Noisy Room

Imagine you are trying to tune a radio to a specific station (the spectral parameter, or $\lambda$). The radio is sitting in a room filled with strange, shifting echoes and obstacles (the potential, or $Q(x)$).

In physics and math, the Dirac equation is the set of rules that describes how a particle (like an electron) moves through this room. Usually, if the room is empty, it's easy to predict where the particle goes. But if the room is full of complex obstacles, predicting the particle's path becomes a nightmare.

For decades, mathematicians have tried to build a "map" to solve this. Previous maps were either:

1. Too slow: You had to calculate the path from scratch for every single radio station you wanted to tune to.
2. Too messy: The map worked well for some stations but got blurry and inaccurate for others.

This paper presents a new, super-efficient map. It says: "We can describe the particle's path using a special stack of building blocks called Neumann series of Bessel functions."

The Core Idea: The "Universal Translator"

The secret weapon in this paper is something called a Transmutation Operator.

Think of the Transmutation Operator as a Universal Translator or a Magic Lens.

• Without the lens: The room is chaotic. The math is hard.
• With the lens: The room looks empty and simple. The math becomes easy.

The authors realized that this "Magic Lens" (which turns the complex problem into a simple one) has a hidden structure. They discovered that this lens is made of a specific type of pattern, much like how a complex song can be broken down into a sequence of simple musical notes.

In this paper, those "notes" are Bessel functions (a specific type of wave pattern used in physics).

The Breakthrough: A "Lego" Approach

The authors didn't just find the lens; they figured out how to build it using Lego bricks.

1. The Bricks (Bessel Functions): These are the standard, well-understood shapes. Everyone knows how they behave.
2. The Instructions (Recursive Integrals): The paper provides a simple, step-by-step recipe to figure out exactly how many of each Lego brick you need and how to stack them. You don't need to guess; you just follow the recipe.
3. The Result (Uniform Convergence): This is the most exciting part. In previous methods, if you wanted to tune to a very high-frequency station (a large $\lambda$), your map would start to crumble and become inaccurate.
   • The Old Way: Like trying to draw a perfect circle with a ruler; it gets messier the bigger the circle gets.
   • This New Way: Like using a compass. No matter how big the circle (or how high the frequency), the accuracy stays perfect. The "error" doesn't get worse as you go further out.

Why Does This Matter? (The "Super-Computer" Effect)

The authors tested their new method on a computer. Here is what happened:

• Speed: They could calculate hundreds of different particle states (eigenvalues) almost instantly.
• Accuracy: Even for the most extreme, high-energy states, the computer got the answer right down to the very last decimal place (machine precision).
• Versatility: Because the method is so clean, it can be used not just to predict where a particle goes (Direct Problem), but also to figure out what the room looks like just by watching the particle (Inverse Problem). This is crucial for things like medical imaging or understanding the structure of materials.

The "Zakharov-Shabat" Connection

The paper also mentions a related system called the Zakharov-Shabat system. Think of this as a "twin" of the Dirac equation. Because the authors built such a robust map for the Dirac equation, they could easily adapt their Lego instructions to map the twin system as well. This opens the door to solving complex problems in fiber optics and laser physics that were previously very difficult to compute.

Summary in a Nutshell

• The Problem: Solving complex particle equations is usually slow and loses accuracy at high energies.
• The Solution: The authors found a way to break the solution down into a stack of simple, predictable waves (Bessel functions).
• The Magic: They created a recipe to build these waves that works perfectly for all energy levels, not just the easy ones.
• The Benefit: Scientists can now simulate complex quantum systems faster and more accurately than ever before, allowing them to solve problems that were previously too computationally expensive.

It's like discovering that instead of manually calculating the trajectory of every single raindrop in a storm, you can just describe the whole storm using a single, perfect mathematical formula that works no matter how hard it rains.

Technical Summary

1. Problem Statement

The paper addresses the numerical and analytical solution of the one-dimensional Dirac equation in canonical form:
$$B \frac{dY}{dx} + Q(x)Y = \lambda Y, \quad x \in (0, b)$$
where $B = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}$, $Q(x)$ is a $2 \times 2$ matrix potential with entries in $L^2([0, b], \mathbb{C})$, and $\lambda$ is a complex spectral parameter.

The authors aim to find a representation of the matrix solution $U(\lambda, x)$ (satisfying $U(\lambda, 0) = I$) that:

1. Is uniformly convergent with respect to the spectral parameter $\lambda$.
2. Allows for the efficient computation of large sets of eigenvalues and eigenvectors without accuracy deterioration for high eigenvalues.
3. Provides an exact representation (rather than a pure approximation) derived from the properties of the transmutation operator.

2. Methodology

The core methodology relies on the theory of transmutation operators and their expansion into Neumann series of Bessel functions (NSBF).

A. Transmutation Operator

The authors utilize a transmutation operator $T$ that maps solutions of the free Dirac equation ($Q=0$) to solutions of the perturbed equation ($Q \neq 0$).
$$U(\lambda, x) = U_0(\lambda, x) + \int_{-x}^{x} K(x, t) U_0(\lambda, t) \, dt$$
where $U_0(\lambda, x)$ is the fundamental solution for $Q=0$, and $K(x, t)$ is the integral kernel satisfying a Goursat problem (a hyperbolic PDE with specific boundary conditions).

B. Fourier-Legendre Expansion

Instead of solving the Goursat problem directly (which is computationally expensive), the authors expand the kernel $K(x, t)$ using Legendre polynomials $P_n$:
$$K(x, t) = \sum_{n=0}^{\infty} \frac{1}{x} K_n(x) P_n\left(\frac{t}{x}\right)$$
This expansion is justified by the smoothness of the kernel $K$ relative to the smoothness of the potential $Q$.

C. Derivation of the NSBF Representation

By substituting the Legendre expansion into the integral representation of $U(\lambda, x)$ and utilizing the orthogonality of Legendre polynomials and known integral identities involving spherical Bessel functions $j_n(z)$, the authors derive the main representation:
$$U(\lambda, x) = U_0(\lambda, x) + 2 \sum_{n=0}^{\infty} (-1)^n K_{2n}(x) j_{2n}(\lambda x) - 2 \sum_{n=0}^{\infty} (-1)^n K_{2n+1}(x) B j_{2n+1}(\lambda x)$$
This is the Neumann Series of Bessel Functions (NSBF).

D. Recursive Computation of Coefficients

A critical contribution is the derivation of a recursive system of differential equations to compute the matrix coefficients $K_n(x)$ (or transformed variables $\theta_n(x) = x^n K_n(x)$) without solving the PDE for $K$.
The coefficients are generated via:

1. Base cases: $K_0(x)$ is derived from the solution of the homogeneous equation ($B Y' + QY = 0$).
2. Recursion: For $n \geq 1$, $\theta_n(x)$ is computed using an integral operator $S$ acting on previous coefficients $\theta_{n-1}$ and $\theta_{n-2}$.
   $$\theta_n(x) = \frac{2n+1}{2n-3} \left[ x^2 \theta_{n-2} + S\left( -(2n-1)x B \theta_{n-2} + (2n-3)\theta_{n-1}B \right) \right]$$
   This allows for the sequential calculation of coefficients using standard ODE solvers.

3. Key Contributions

1. Exact NSBF Representation for Dirac Systems: The paper extends the NSBF method, previously applied to Schrödinger equations, to the $2 \times 2$ Dirac system. This provides a unified framework for solving direct and inverse spectral problems.
2. Uniform Convergence: The authors prove that the truncated series converges uniformly with respect to $\lambda$. Specifically, the error bound depends on the smoothness of $Q$ but not on $\lambda$ (for real $\lambda$), allowing for the computation of arbitrarily high eigenvalues with constant accuracy.
3. Recursive Algorithm: They provide an explicit, efficient recursive algorithm to compute the series coefficients. This avoids the need to solve initial value problems for every spectral parameter $\lambda$, a major bottleneck in traditional numerical methods (like shooting methods).
4. Extension to ZS-AKNS Systems: The method is generalized to the Zakharov-Shabat (ZS-AKNS) system, which is fundamental in the inverse scattering method for non-linear Schrödinger equations.
5. Connection to Formal Powers: The appendix establishes a rigorous link between the NSBF coefficients and the "formal powers" used in the Spectral Parameter Power Series (SPPS) method, unifying two distinct analytical approaches.

4. Results and Numerical Performance

• Convergence Analysis: The paper establishes error bounds for the truncated series $U^N(\lambda, x)$. If $Q \in C^{p+1}$, the error decays as $O(N^{-(p+1)})$. If $Q$ is in a Sobolev space $W^{p+1, 2}$, similar decay rates are proven using $L^2$ norms.
• Numerical Example: The authors tested the method on a spectral problem with boundary conditions $u(0)=u(1)=0$.
  • They computed 240 eigenvalues (ranging from -105 to 134) using only $N=16$ terms in the series.
  • Accuracy: For large eigenvalues, the computed values coincided with exact values up to machine precision ($10^{-14}$).
  • Efficiency: The method demonstrated that accuracy does not deteriorate as the spectral parameter increases, a common failure point for other numerical methods.

5. Significance

• Computational Efficiency: The method drastically reduces computational cost for spectral problems. Once the coefficients $K_n(x)$ are computed (which depends only on $Q$), the solution $U(\lambda, x)$ can be evaluated for any $\lambda$ instantly. This is superior to methods requiring re-integration for each $\lambda$.
• Inverse Problems: The representation opens the door to solving inverse spectral problems for the Dirac equation numerically, a task that has been challenging due to the lack of stable, high-accuracy forward solvers for large spectral ranges.
• Robustness: The uniform convergence property makes the method particularly suitable for problems involving high-energy states or large spectral parameters, where traditional discretization or shooting methods often fail or require excessive computational resources.

In summary, this paper presents a powerful analytical and numerical tool that transforms the solution of the 1D Dirac equation into a rapidly convergent series, enabling high-precision computation of spectral data across the entire complex plane.