Investigating a two-level algorithm for fermionic observables

We investigate the combination of a two-level sampling algorithm with distillation techniques to compute disconnected fermionic correlation functions. The method relies on a factorization of the quark propagator into domain-local contributions that depend only on the gauge fields within overlapping temporal regions, enabling independent submeasurements of each term through a two-level sampling strategy. The statistical uncertainty follows the expected inverse scaling with the number of second-level submeasurements, up to exponential boundary effects, in contrast to the inverse square-root scaling of standard sampling. This leads to an exponential reduction of statistical errors at negligible additional computational cost compared to the gain in statistical precision. The method is tested on pure gauge ensembles, providing a controlled benchmark for its forthcoming application to dynamical QCD studies.