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 two-level estimators exhibit the expected $1/N_1^2$ scaling of the variance, up to exponential boundary effects, and achieve an exponential reduction of statistical errors at nearly the same computational cost as standard sampling. The method is tested on pure gauge ensembles, providing a controlled benchmark for its forthcoming application to dynamical QCD studies of glueball and isosinglet meson correlation functions.