Regularized Unfolding of gamma-ray Spectra for Nuclear Physics Applications

TL;DR

This work presents a Regularized Maximum Likelihood Estimation (RMLE) framework for unfolding gamma-ray spectra within the Oslo Method, treating unfolding as a nonnegative Poisson inverse problem with explicit modeling of background and contaminants. By embedding the detector response in a forward model and applying physically constrained regularization, RMLE yields smoother, more robust spectra with calibrated uncertainty intervals, and can be computed efficiently on GPUs. The authors develop a comprehensive treatment of ill-posedness, null-space degeneracy, and uncertainty quantification via Monte Carlo resampling, and they demonstrate favorable performance for low-complexity spectra while outlining the limits for high-complexity cases. Compared to the traditional Folding Iteration with Compton Subtraction (FICS), RMLE provides principled uncertainty estimates and reduced overfitting, offering a practically impactful path for reliable Oslo-Method analysis and broader inverse problems in nuclear spectroscopy.

Abstract

Reconstructing gamma-ray spectra from detector measurements is an ill-posed inverse problem. Standard methods, such as Folding Iteration with Compton Subtraction (FICS), provide point estimates but lack calibrated uncertainties and may bias the spectrum. We introduce an unfolding framework based on regularized maximum-likelihood estimation (RMLE) that enforces non-negativity and detector-response constraints while explicitly modeling background and contaminant contributions. Simulations and analytical results show that RMLE yields smoother reconstructions with well-calibrated confidence intervals and outperforms existing techniques for low-complexity spectra. Although high-complexity data remain challenging, the intervals produced by RMLE maintain correct coverage.

Figures (36)

• Figure 1: (a) Each component of the discrete response for true $E_{\gamma}=9MeV$. The peaks have been scaled down to make them visually comparable. (b) A single sharp peak at 9MeV folded with $\mathbf{G}_{\gamma}\mathbf{D}$.
• Figure 2: The discrete OSCAR response $\mathbf{D}$ for a $1000\times 1000$ matrix. The discrete peaks $\boldsymbol{p}_{\text{f}}, \boldsymbol{p}_{\text{s}}$ and $\boldsymbol{p}_{\text{d}}$ are along the diagonal and at offsets 511keV and 1024keV, respectively. The $\boldsymbol{p}_{\text{a}}$ is a sharp vertical structure at measured $E_{\gamma} = 511keV$. The remaining bulk is the $\mathbf{P}_{\text{c}}$ component. The inset axes shows an example of a single row at true $E_{\gamma} = 5MeV$. The color scale is logarithmic, scaled to prevent outlying bins from affecting the color. The numbers on the colorbar indicate $0.3\%$ of the bins lie above the range.
• Figure 3: Condition numbers versus matrix order for matrices $\mathbf{D}$, $\mathbf{G}_{\gamma}{}$, $\mathbf{G}_{\gamma}{}\mathbf{D}$, $\mathbf{G_{\text{in}}}$ and $\mathbf{G}_{\gamma}\mathbf{D}\mathbf{G_{\text{in}}}$. The smoothing operators $\mathbf{G}_{\gamma}{}$ and $\mathbf{G_{\text{in}}}$ exhibit substantially higher condition number growth with increasing order than $\mathbf{D}$. Shaded bands represent the range of condition numbers across 100 instances with small perturbations, simulating numerical fluctuations in matrix construction. The resolution is held constant at $\sigma_{\text{in}} = 40keV$ and $\sigma_\gamma(1330keV) = 40keV$. The energy range is from $0$ to 10MeV, with $\Delta E$ determined by the order.
• Figure 4: Examples of smoothing matrices for the initial-excitation-energy axis, $\mathbf{G_{\text{in}}}$ (a), and gamma-ray-energy axis, $\mathbf{G}_{\gamma}$ (b). The inset axes show examples of the Gaussians at specific true energies marked with dotted lines. The color scale (not shown) is logarithmic with $0$ mapped to white.
• Figure 5: (a) Condition number and (b) mean row distance of Gaussian smoothing matrix $\mathbf{G}$ versus resolution parameter $\sigma$, with fixed bin width $\Delta E = 10keV$. The dashed line in the bottom panel marks $\sqrt{2}$. The shaded bands represent the range of condition numbers across 100 instances with small perturbations, simulating numerical fluctuations in matrix construction. The transition from constant to linear distance scaling occurs when $\sigma \approx \Delta E$, coinciding with rapid growth in condition number. For $\mathbf{G}_{\gamma}$, the $\sigma_\gamma(E_{\gamma})$ is calibrated so that the mean $\sigma_\gamma$ over $E_{\gamma}$ equals $\sigma$.

Theorems & Definitions (11)

• Definition 2.1: Poisson Distribution
• Definition 2.2: Maximum Likelihood
• Definition 2.3: Point Measure
• Definition 2.4: Hadamard Criteria
• Theorem 2.1
• Corollary 2.2
• Definition 2.5
• Proposition 2.3
• Theorem 2.4
• Proposition A.1