A general recursion for integrals involving products of Hermite polynomials and its applications

1. Problem Statement

The paper addresses a critical computational challenge in the ab initio simulation of correlated few-body quantum systems confined in one-dimensional (1D) harmonic traps. These simulations, often utilizing the Configuration Interaction (CI) approach, require the evaluation of matrix elements for two-body and three-body contact interactions.

These matrix elements involve integrals of products of $N$ Hermite polynomials ($H_n(x)$) weighted by a Gaussian function. Specifically, the authors focus on three integral types:

• $W_{ijk\ell}$: Two-body interaction matrix elements (product of 4 Hermite polynomials, $N=4$).
• $Y_{ijk\ell}$: Momentum-dependent two-body interaction matrix elements (involving derivatives of Hermite polynomials).
• $U_{ijk\ell mn}$: Three-body interaction matrix elements (product of 6 Hermite polynomials, $N=6$).

Existing Limitations:

• Numerical Instability: Traditional analytical formulas for these integrals rely heavily on explicit factorials of large integers. As the polynomial degree ($M$) increases to ensure numerical convergence, these factorials cause numerical overflows and instabilities.
• Computational Inefficiency: Symbolic math toolboxes become prohibitively slow for large $M$ due to the polynomial growth in the number of distinct integrals.
• Quadrature Sensitivity: Direct numerical integration (quadrature) is highly sensitive to grid point selection and integration intervals, often failing to achieve high precision for high-degree polynomials.
• Lack of Generalization: While some analytical solutions exist for $W_{ijk\ell}$, they are algebraically complex and difficult to implement efficiently. No analytical expressions existed in the literature for $Y_{ijk\ell}$ or $U_{ijk\ell mn}$.

2. Methodology

The authors derive a general recursive formula for the integral of the product of $N$ Hermite polynomials, denoted as $T_{a_1 \dots a_N}$. The derivation relies on two fundamental mathematical tools:

1. Recurrence Relations of Hermite Polynomials: Utilizing the relation $H_n(x) = 2xH_{n-1}(x) - 2(n-1)H_{n-2}(x)$.
2. Integration by Parts: Applying this technique to the Gaussian-weighted integral to reduce the order of the polynomials.

Key Derivation Steps:

• The integral $T_{a_1 \dots a_N}$ is rewritten by substituting the recurrence relation for one of the polynomials.
• Integration by parts is applied to the term containing $x$, utilizing the derivative property $\frac{d}{dx}H_n(x) = 2nH_{n-1}(x)$ and the derivative of the Gaussian weight.
• This process yields a recursive relation where an integral of order $V = \sum a_i$ is expressed as a linear combination of integrals of order $V-2$.

The Recursive Formula:
$$T_{a_1 \dots a_N} = \left( \frac{2}{N} - 1 \right) \sqrt{\frac{a_1 - 1}{a_1}} T_{(a_1-2)a_2 \dots a_N} + \frac{2}{N} \sum_{i=2}^{N} \sqrt{\frac{a_i}{a_1}} T_{(a_1-1)\dots(a_i-1)\dots a_N}$$

Selection Rule:
Due to the symmetry of the integration domain $(-\infty, +\infty)$ and the even nature of the Gaussian weight, the integral is non-zero if and only if the sum of the indices $V = \sum_{i=1}^N a_i$ is an even integer. If $V$ is odd, the integral vanishes. This allows the algorithm to skip zero-valued calculations efficiently.

Base Case:
The recursion terminates at the lowest level ($V=0$), where all indices are zero. This reduces to a fundamental Gaussian integral with a known closed-form analytical solution, eliminating the need for any numerical quadrature.

3. Key Contributions

• Factorial-Free Formulation: The most significant contribution is the complete elimination of explicit factorials from the calculation. By using square roots of ratios of indices, the method avoids the numerical overflow and instability associated with large integer factorials.
• Generalization: The derivation provides a unified framework for any number of polynomials ($N \geq 2$), covering the specific cases of $N=4$ ($W$), $N=6$ ($U$), and the derivative-based $Y$.
• Analytical Expressions for New Integrals: The paper provides the first recursive formulas for $Y_{ijk\ell}$ and $U_{ijk\ell mn}$, which were previously unavailable in the literature.
• Dynamic Programming Implementation: The recursive structure is optimized for computer implementation using dynamic programming (caching intermediate results), ensuring high efficiency.
• Open-Source Tools: The authors provide a Mathematica notebook and a Python script implementing these recursions, available via GitHub.

4. Results

The authors validated their method through extensive numerical testing:

• Precision: The Python implementation (using double-precision floating-point arithmetic) was compared against Mathematica's symbolic integration (using 500 significant digits).
• Accuracy: The results showed perfect agreement at the level of machine precision (absolute errors on the order of $10^{-17}$), even for extremely high-order terms (e.g., $W_{500,500,500,500}$).
• Scalability: The method successfully computed integrals for indices up to $M=500$ (and higher), a regime where symbolic integration becomes impractically slow and factorial-based methods fail due to overflow.
• Efficiency: The recursive approach is significantly faster than brute-force quadrature or symbolic expansion, as it relies solely on elementary arithmetic operations.

5. Significance and Applications

• Quantum Physics: The method enables high-precision ab initio simulations of few-body systems in 1D harmonic traps. This is crucial for studying phenomena like fermionization, anyonic systems, and quantum chaos in ultracold atomic gases.
• Computational Physics: It offers a robust, stable, and efficient alternative to existing methods for calculating matrix elements in Configuration Interaction approaches.
• Broader Mathematical Impact: While focused on Hermite polynomials, the authors note that the derivation technique can be extended to other classes of classical orthogonal polynomials.
• Practical Utility: By providing a stable, factorial-free algorithm, this work removes a major bottleneck in the simulation of strongly correlated quantum systems, allowing researchers to push the boundaries of system size and interaction complexity.

In summary, this paper presents a mathematically elegant and computationally robust solution to a long-standing numerical problem in quantum many-body physics, replacing unstable analytical formulas and inefficient numerical integration with a stable, recursive, and highly accurate algorithm.