Fourier Transforms of Color Glass Condensate Multi-Wilson-Line Correlators via Filon Quadrature

Abstract: Calculating cross sections in the Color Glass Condensate effective theory requires Fourier transforms of multi-Wilson-line correlators from transverse coordinate space to transverse momentum space. Under the common assumption of impact-parameter independence, each transform reduces to a set of Hankel transforms whose Bessel-function kernels oscillate rapidly at phenomenologically relevant momenta, making direct quadrature prohibitively expensive. We present a Filon-type quadrature, applicable to any integrand, that integrates these oscillatory factors in closed form on the stored coordinate grid, reducing each Hankel transform to a precomputed weight vector and the full nested transform chain to a sequence of matrix products. We develop and validate the method on the deep inelastic scattering dijet cross section beyond the correlation-limit approximation, where an exprel-based reformulation of the quadrupole Wilson-line correlator removes a numerical $0/0$ instability inherent to its standard parametrization. Porting the calculation to the Graphics Processing Unit (GPU), with custom CUDA kernels that fuse the momentum-space contraction directly into the correlator evaluation, brings the runtime for one dipole input down to about two minutes on a single NVIDIA A800, from several hours on a multi-core Central Processing Unit (CPU). We further generalize the algorithm to three sequential Hankel transforms and validate the resulting six-dimensional transform against an analytic Gaussian integrand family with closed-form results at every stage. This general, process-independent algorithm is directly applicable to next-to-leading-order proton-nucleus and electron-ion scattering cross-section calculations performed without the correlation-limit approximation. The code is publicly available at https://github.com/CCNU-CGC-py/FFT_filon.

• The paper develops a Filon-quadrature method that converts nested oscillatory Hankel transforms into matrix-weight contractions, enabling efficient multi-Wilson-line Fourier transforms beyond the correlation-limit approximation.
• The GPU implementation reduces leading-order DIS dijet calculations from hours on multicore CPUs to 49–84 seconds on one A800 while preserving CPU results through deep cancellation minima and high-resolution angular convergence tests.
• The method generalizes to three sequential Hankel transforms for six-dimensional Fourier integrals, achieving better than 4×10⁻⁶ accuracy on analytic Gaussian tests and supporting future NLO CGC phenomenology.

Motivation and problem statement

Cross sections in the Color Glass Condensate (CGC) effective theory are formulated in transverse coordinate space, where small-$x$ evolution of Wilson-line correlators is simplest, but are measured differentially in transverse momentum. Their evaluation therefore requires nested two-dimensional Fourier transforms of multi-Wilson-line correlators. Under the standard assumption of impact-parameter independence, each transform reduces to a Hankel (Bessel) transform; at phenomenologically relevant momenta the Bessel kernels oscillate with $q\,b$ reaching $\mathcal{O}(10^3)$ across the coordinate grid, so direct quadrature becomes prohibitively expensive. Prior numerical treatments of the four-point (quadrupole) correlator were limited to relatively low transverse momenta (Lappi et al., 2012, Mäntysaari et al., 2019), while the Factorized Dipole Approximation bypasses the transform entirely but assumes uncorrelated color fluctuations and is valid only in the regime $Q_s^2 S_\perp \gg 1$, discarding the connected multi-gluon correlations relevant to saturation.

The paper develops a Filon-type quadrature [filon1930iii] that integrates the oscillatory Bessel factors in closed form on a stored coordinate grid, reducing each Hankel transform to a precomputed weight vector and the full nested transform chain to matrix products. The method is developed and validated on the leading-order DIS dijet cross section beyond the correlation-limit approximation, then generalized to three sequential Hankel transforms — a genuinely six-dimensional Fourier transform — validated against an analytic Gaussian family.

The dijet observable and the exprel reformulation

The observable is inclusive dijet production in DIS, $\gamma^*(q_\gamma) + p/A \to q\bar q + X$, in the dipole picture. The cross section involves the quadrupole correlator $S^{(4)}(\mathbf r,\mathbf r',\mathbf b)$ built from dipole $S$-matrices via the Gaussian-approximation, large-$N_c$ reduction, requiring only the two-point function $S(r)$ as input — either GBW or a tabulated collinearly-improved BK solution extracted with a neural network (Dai et al., 9 Mar 2026, Dai et al., 30 Jul 2026). The integral carries two Fourier phases, conjugate to the tagged-parton momentum $\mathbf k_1$ and the pair recoil $q\,b$0.

A technical contribution worth noting: the standard parametrization of $q\,b$1 takes the form $q\,b$2 on the codimension-one surface where $q\,b$3. Rewriting the numerator as $q\,b$4 yields $q\,b$5 with $q\,b$6, which removes this instability exactly, using a Taylor guard for $q\,b$7.

Angular reduction and Filon weight vectors

Expanding both plane-wave phases in Bessel series and integrating over the global angle enforces the selection rule $q\,b$8, collapsing the six-fold integral into a sum over azimuthal harmonics of double Hankel transforms $q\,b$9, truncated at $\mathcal{O}(10^3)$0 once the mode sum saturates.

The Filon construction proceeds interval by interval: the smooth factor $\mathcal{O}(10^3)$1 is fit by a degree-4 polynomial through five grid nodes via Lagrange cardinal polynomials (the Vandermonde system is well-conditioned after rescaling), and the oscillatory factor is integrated against each monomial in closed form using Taylor-expanded local moments,

$\mathcal{O}(10^3)$2

By Stirling's estimate the series converges beyond $\mathcal{O}(10^3)$3; a common truncation order $\mathcal{O}(10^3)$4 (capped at 400), set by the widest interval, guarantees convergence everywhere. Regrouping the sums expresses each transform as $\mathcal{O}(10^3)$5: the integrand enters only linearly through its sampled values, and the weight vector depends solely on the (grid, momentum, order) triple — reusable across polarizations and dipoles. Stacking rows over momenta turns each transform into a matrix product.

CPU reference implementation

The CPU implementation separates into stage 1, building the momentum-free inner table $\mathcal{O}(10^3)$6 of shape $\mathcal{O}(10^3)$7, and stage 2, applying the two Filon contractions. With $\mathcal{O}(10^3)$8, $\mathcal{O}(10^3)$9, $Q_s^2 S_\perp \gg 1$0, stage 1 dominates: about 5.5 hours on 16 cores at $Q_s^2 S_\perp \gg 1$1 GeV and roughly 22 hours (scaled) at $Q_s^2 S_\perp \gg 1$2 GeV, versus tens of seconds for stage 2. Convergence is verified in three ways: saturation of the angular mode sum, grid convergence of the Filon contraction, and visual inspection of tables and cross-section slices. A limitation is identified here: at $Q_s^2 S_\perp \gg 1$3 GeV the default angular resolution leaves some production points unconverged, even producing negative truncated reconstructions; resolving them at $Q_s^2 S_\perp \gg 1$4 would multiply the dominant $Q_s^2 S_\perp \gg 1$5 scaling by sixteen, making a CPU study impractical — a key motivation for the GPU port.

GPU implementations

Two CUDA implementations preserve the mathematics exactly. The table builder ports stage 1 to GPU with two optimizations: transcendental functions ($Q_s^2 S_\perp \gg 1$6, regularized $Q_s^2 S_\perp \gg 1$7, $Q_s^2 S_\perp \gg 1$8) are precomputed as spline-coefficient tables evaluated by polynomial interpolation, and each integrand factor is computed only on its actual variable dependence via four custom kernels plus a cuBLAS angular projection. Nsight Compute profiling shows all kernels compute-bound (SM 83–95%, DRAM 5–11%), confirming the Roofline interpretation that FP64 throughput, not memory bandwidth, limits the calculation. On an A800 with full FP64 rate, stage 1 drops from 167 s to 23 s at $Q_s^2 S_\perp \gg 1$9 GeV (speedup ≈ 7.3 over the RTX 3090) and from 666 s to 74 s at $\gamma^*(q_\gamma) + p/A \to q\bar q + X$0 GeV (≈ 9.0). GPU cross sections overlay the CPU reference, including through deep cancellation minima.

The high-resolution stress test resolves the earlier limitation conclusively: at the previously unconverged point ($\gamma^*(q_\gamma) + p/A \to q\bar q + X$1 GeV, $\gamma^*(q_\gamma) + p/A \to q\bar q + X$2 GeV, $\gamma^*(q_\gamma) + p/A \to q\bar q + X$3 GeV), raising resolution to $\gamma^*(q_\gamma) + p/A \to q\bar q + X$4, $\gamma^*(q_\gamma) + p/A \to q\bar q + X$5 recovers a positive angular minimum, with the longitudinal partial sum at $\gamma^*(q_\gamma) + p/A \to q\bar q + X$6 differing from $\gamma^*(q_\gamma) + p/A \to q\bar q + X$7 by only $\gamma^*(q_\gamma) + p/A \to q\bar q + X$8 relative. The apparent failure at default resolution is therefore purely an angular-truncation effect, not an algorithmic one.

The fused implementation removes the stage boundary altogether: since both Hankel transforms are linear sums over table entries, the $\gamma^*(q_\gamma) + p/A \to q\bar q + X$9-contraction is applied incrementally as $S^{(4)}(\mathbf r,\mathbf r',\mathbf b)$0 values are computed, and the full inner table never exists. Peak memory drops to ~10 MB intermediates plus a few MB output (~$S^{(4)}(\mathbf r,\mathbf r',\mathbf b)$1 rows instead of multi-gigabyte tables). End-to-end times reach 49 s at $S^{(4)}(\mathbf r,\mathbf r',\mathbf b)$2 GeV and 84 s at $S^{(4)}(\mathbf r,\mathbf r',\mathbf b)$3 GeV on one A800 — down from several hours on a multi-core CPU. Agreement with the table-builder path is at the level of absolute differences $S^{(4)}(\mathbf r,\mathbf r',\mathbf b)$4–$S^{(4)}(\mathbf r,\mathbf r',\mathbf b)$5 and relative differences below $S^{(4)}(\mathbf r,\mathbf r',\mathbf b)$6 over all modes, momenta, and polarizations, attributable to floating-point reordering. One trade-off is explicit: because Filon weights encode the momenta, changing $S^{(4)}(\mathbf r,\mathbf r',\mathbf b)$7 requires recomputation, forfeiting the stored table's reusability at new momenta.

Generalization to six dimensions

The general case attaches independent momenta to all three coordinate variables, giving a six-dimensional transform evaluated as three sequential Hankel transforms with selection rule $S^{(4)}(\mathbf r,\mathbf r',\mathbf b)$8, leaving two signed mode indices. Because the momentum-free intermediate would be orders of magnitude too large to store per Bjorken-$S^{(4)}(\mathbf r,\mathbf r',\mathbf b)$9 and $S$0 combination, only the fused streaming form is viable; if the intermediate exceeds memory, the calculation is chunked over groups of $S$1 values. Validation uses an analytic positive-definite Gaussian family whose transform is known in closed form: numerical and exact results agree to better than $S$2 at every tested point, verifying the angular projection, all three Hankel transforms, and the final mode reconstruction.

Limitations and open questions

The paper is explicit about what remains undone. The three-transform algorithm has not been applied to any physical observable; its performance study with the analytic test function is not representative of process-specific calculations, since each process requires its own CUDA kernels, memory layouts, grids, and optimization. Production deployment is deferred to next-to-leading-order proton-nucleus and electron-ion cross sections without the correlation-limit approximation. Additionally, the kernel decomposition exploiting the factor structure of the DIS integrand is specific to that observable, and the fused path's momentum rigidity means new momentum grids require full recomputation.

Conclusion

This work provides a general, process-independent algorithm for the nested Hankel transforms central to CGC phenomenology: closed-form Filon moments collapse oscillatory Bessel integrals into reusable weight vectors, an exprel identity eliminates the quadrupole's $S$3 instability, and staged validation — CPU reference, CUDA table builder, fused streaming implementation, analytic six-dimensional benchmark — establishes correctness at each step. The practical outcome is a reduction of the DIS dijet computation from hours on multi-core CPUs to one to two minutes on a single commercial GPU, with the code publicly available, providing a direct foundation for forthcoming NLO $S$4 and electron-ion applications beyond the correlation limit.