Free-particle Green's function matrix elements over spherical Gaussian and plane-wave-modulated Gaussian basis functions

Abstract: Free-particle Green's function plays a central role in the theoretical description of electron scattering and autoionization processes in quantum physics and chemistry. Recently, Gaussian basis set approaches have become increasingly important in applications to unbound and metastable electronic states. However, the practical use of such methods has been limited by the lack of efficient and compact analytical expressions for matrix elements of the free-particle Green's function in Gaussian-based representations. Here we present a novel, general analytical framework for the evaluation of one- and two-center matrix elements of the free-particle Green's function over spherical Gaussian basis functions and plane-wave-modulated spherical Gaussians. The derivation is based on the Fourier transform of Gaussian functions together with the addition theorem of harmonic polynomials, leading to compact closed-form expressions and efficient recurrence relations. We also analyze the asymptotic behavior of the free-particle Green's function matrix elements, which is essential in the description of low-energy continuum electrons using finite Gaussian basis sets.

• The paper develops compact analytical expressions and angular-momentum recurrences for outgoing free-particle Green’s function matrix elements over spherical Gaussian basis functions, reducing all two-center cases to reusable radial seed integrals.
• The method extends to plane-wave-modulated spherical Gaussians through complex Gaussian centers, while asymptotic complementary-error-function expansions stabilize calculations for very diffuse functions.
• The implementation in Libqints and Q-Chem agrees with high-precision references and prior Cartesian-Gaussian results, and exponent scans suggest matching Gaussian widths to the electron de Broglie wavelength when designing continuum bases.

Motivation and context

Processes involving free electrons—electron–molecule scattering, photoionization, Auger decay, and Penning ionization—require a rigorous treatment of continuum electronic states that are not square-integrable and therefore cannot be represented exactly in finite $L^2$ basis sets. Gaussian-type orbitals (GTOs) dominate bound-state quantum chemistry because of the efficient integral technology originating from Boys' work, but their use for continuum problems has been hampered by one specific technical obstacle: the absence of compact, general analytical expressions for matrix elements of the free-particle Green's operator

$$\hat{G}_0^{(+)}(E) = \lim_{\epsilon \to 0^+} \frac{1}{E - \hat{T} + i\epsilon},$$

which is central to the Lippmann–Schwinger equation. Prior Cartesian-Gaussian formulations by Östlund (2605.18564), Levin et al., and Čársky et al. were formally exact but algebraically cumbersome at high angular momenta and did not reuse intermediate quantities efficiently across symmetry-related integrals, which the authors identify as the reason these methods saw little practical adoption.

General momentum-space framework

The paper derives matrix elements over spherical GTOs (SGTOs), $$\phi_{lm}(\alpha,\mathbf{r}) = N_l(\alpha)e^{-\alpha r^2}\mathcal{Y}_{lm}(\mathbf{r})$$, by inserting momentum-space completeness relations and exploiting the fact that the Green's function is diagonal in momentum space:

$$\langle \mathbf{p}|\hat{G}^{(+)}(k)|\mathbf{q}\rangle = \lim_{\epsilon \to 0^+} \frac{\delta(\mathbf{p}-\mathbf{q})}{k_0^2 - q^2 + i\epsilon}.$$

After expanding the product of solid harmonics via Gaunt coefficients and using the partial-wave expansion of the plane-wave factor, every two-center matrix element reduces to an angular sum weighted by radial integrals of the form

$$I_l^{P'}(\eta,R,k_0) = \int_0^\infty \frac{q^{P'} e^{-\eta q^2} j_l(qR)}{k_0^2 - q^2 + i\epsilon}\,dq,$$

with $P' = 2 + l_a + l_b$. The same structure holds for real spherical harmonics with real Gaunt coefficients, making the framework directly compatible with production quantum chemistry codes. One-center integrals ($R=0$) collapse to a diagonal form $\delta_{l_al_b}\delta_{m_am_b}$.

Recurrence relations

The computational core of the work is a set of stable recurrences. A Bessel-function recurrence,

$$I_{l+1}^{P'+1} = \frac{2l+1}{R}I_l^{P'} - I_{l-1}^{P'+1},$$

means only $l=0$ and $$\hat{G}_0^{(+)}(E) = \lim_{\epsilon \to 0^+} \frac{1}{E - \hat{T} + i\epsilon},$$0 seed integrals must be computed; all higher angular momenta follow recursively. The seeds themselves are generated from closed-form expressions involving the complementary error function, combined with differentiation with respect to $$\hat{G}_0^{(+)}(E) = \lim_{\epsilon \to 0^+} \frac{1}{E - \hat{T} + i\epsilon},$$1. For the one-center case, the fundamental integral is expressed through the Faddeeva function, $$\hat{G}_0^{(+)}(E) = \lim_{\epsilon \to 0^+} \frac{1}{E - \hat{T} + i\epsilon},$$2, with higher powers obtained by a binomial recurrence. Analogous kinetic-energy recurrences (Appendix) serve as validation targets since both operators share identical angular structure.

Plane-wave-modulated SGTOs

The formalism extends naturally to plane-wave-modulated Gaussians (PW-SGTOs), $$\hat{G}_0^{(+)}(E) = \lim_{\epsilon \to 0^+} \frac{1}{E - \hat{T} + i\epsilon},$$3, which embed oscillatory behavior directly into the basis and can drastically reduce the number of functions needed to represent continuum states compared to pure Gaussians. The key observation is that the plane wave shifts the Gaussian center into the complex plane: the modulated $$\hat{G}_0^{(+)}(E) = \lim_{\epsilon \to 0^+} \frac{1}{E - \hat{T} + i\epsilon},$$4-type Gaussian equals a renormalized Gaussian centered at $$\hat{G}_0^{(+)}(E) = \lim_{\epsilon \to 0^+} \frac{1}{E - \hat{T} + i\epsilon},$$5. Applying the addition theorem of solid harmonics then expresses PW-SGTO matrix elements entirely in terms of standard SGTO matrix elements evaluated at the complex displacement $$\hat{G}_0^{(+)}(E) = \lim_{\epsilon \to 0^+} \frac{1}{E - \hat{T} + i\epsilon},$$6, multiplied by angular-coupling coefficients. The real-harmonic version requires four case distinctions depending on the signs of the magnetic indices, each reducing to a nine-term sum. Notably, the authors point out that the earlier prototype PW-Cartesian-Gaussian formulas of Colle et al. contained inconsistencies in overlap normalization factors—an explicit claim of error in prior literature.

Numerical stability and asymptotics

Direct evaluation of the Erfc-based expressions fails in the regime $$\hat{G}_0^{(+)}(E) = \lim_{\epsilon \to 0^+} \frac{1}{E - \hat{T} + i\epsilon},$$7 (very diffuse Gaussians), where cancellation between exponentially large and small terms causes overflow. The paper supplies asymptotic expansions of the complementary error function for this regime, applied above a conservative switching condition $$\hat{G}_0^{(+)}(E) = \lim_{\epsilon \to 0^+} \frac{1}{E - \hat{T} + i\epsilon},$$8 with terms up to $$\hat{G}_0^{(+)}(E) = \lim_{\epsilon \to 0^+} \frac{1}{E - \hat{T} + i\epsilon},$$9. In this regime the matrix elements become purely real for real SGTOs, and the imaginary part vanishes identically—a physically meaningful signature of the loss of on-shell continuum coupling for extremely diffuse basis functions. The large-separation regime $$\phi_{lm}(\alpha,\mathbf{r}) = N_l(\alpha)e^{-\alpha r^2}\mathcal{Y}_{lm}(\mathbf{r})$$0 is shown to be numerically benign, reducing to a simple outgoing-wave form proportional to $$\phi_{lm}(\alpha,\mathbf{r}) = N_l(\alpha)e^{-\alpha r^2}\mathcal{Y}_{lm}(\mathbf{r})$$1.

Validation and numerical behavior

The implementation was developed first in Mathematica for high-precision reference values and then in the Libqints library within Q-Chem. Three independent checks were performed: agreement with Mathematica references, consistency of the shared angular machinery against standard Libqints kinetic-energy integrals, and benchmarking against Čársky et al.'s Cartesian-Gaussian results after transformation to the spherical representation. The paper publishes benchmark tables of one-center, two-center, and PW-SGTO matrix elements to enable independent verification.

The exponent-dependence analysis carries practical weight: the imaginary part of the matrix element exhibits a maximum when the Gaussian width matches the electron de Broglie wavelength, and these maxima shift toward tighter exponents as the electron energy increases (illustrated at 1 eV versus 100 eV). This provides an a priori criterion for constructing even-tempered continuum basis sets concentrated where the Green's function matrix elements vary most strongly.

Limitations and open questions

The paper restricts attention to the outgoing Green's operator, obtaining $$\phi_{lm}(\alpha,\mathbf{r}) = N_l(\alpha)e^{-\alpha r^2}\mathcal{Y}_{lm}(\mathbf{r})$$2 only by Hermitian conjugation, and treats exclusively momentum-diagonal one-electron operators; extension to non-diagonal operators is not addressed. The asymptotic switching threshold $$\phi_{lm}(\alpha,\mathbf{r}) = N_l(\alpha)e^{-\alpha r^2}\mathcal{Y}_{lm}(\mathbf{r})$$3 is described as conservative rather than optimized, and the truncation order $$\phi_{lm}(\alpha,\mathbf{r}) = N_l(\alpha)e^{-\alpha r^2}\mathcal{Y}_{lm}(\mathbf{r})$$4 of the asymptotic series is fixed without a systematic error analysis. For PW-SGTOs, the authors note that exponentially small normalization factors must be folded into the radial integrals to avoid overflow, but a fully robust treatment across all parameter regimes is presented as a practical prescription rather than a proven-stable algorithm. Finally, no end-to-end scattering or autoionization calculation is performed here; applications to Lippmann–Schwinger solvers and ab initio cross sections are deferred to forthcoming publications, so the actual accuracy gains of PW-SGTO representations relative to large pure-Gaussian bases remain quantitatively untested in this work.

Conclusion

This paper provides a complete, angular-momentum-general analytical and recursive framework for free-particle Green's function matrix elements over SGTOs and PW-SGTOs, replacing the cumbersome Cartesian formulations that previously limited practical use. The combination of compact recurrences, asymptotically stabilized evaluation, published benchmark values, and an implementation inside a production integral library positions the formalism for direct incorporation into bound-state electronic structure codes extended to continuum problems. The decisive test—demonstrated accuracy and efficiency in full scattering or resonance calculations—remains open.