Differentiable ABCD matrix solver for Gaussian beamlet decomposition

Gaussian beamlet decomposition is a powerful framework for modeling coherent light propagation in complex optical systems. A key challenge is accurately evaluating the associated ABCD matrix, which is traditionally computed via finite-difference ray tracing with multiple perturbed rays, resulting in numerical errors and sensitivity to perturbation step size. Here, we present a differentiable ABCD matrix solver that enables direct computation of the ABCD matrix without finite-difference approximations. By formulating ray propagation as a differentiable computational graph, the ABCD matrix is obtained via automatic differentiation. We validate the proposed method in both a high-NA microscopic imaging system and a non-imaging beam-shaping system, demonstrating high accuracy and numerical stability compared to conventional approaches. Furthermore, the fully differentiable framework enables gradient-based backpropagation, enabling seamless integration with modern optimization frameworks and inverse wave-optics design.