Modern cosmology relies on the analysis of large galaxy surveys to understand the large-scale structure of the universe. A crucial aspect is the modeling of the power spectrum of biased tracers (galaxies), which do not perfectly reflect the distribution of dark matter. For decades, perturbative templates have been developed in Eulerian and Lagrangian frameworks for the standard cosmological model $\Lambda$CDM. However, to go beyond $\Lambda$CDM and explore modified gravity theories, more sophisticated tools are required that can accurately handle nonlinear regimes.

This work addresses the implementation of the biased perturbative expansion, within the local Lagrangian bias scheme, in the framework of the Hybrid Effective Field Theory (HEFT). HEFT combines perturbation theory with dark matter simulations to model the nonlinear regime. The researchers focused on $f(R)$ gravity, a modified gravity theory that exhibits scale-dependent growth and a "chameleon screening" effect, making it a particularly challenging scenario for Lagrangian perturbation theory calculations and for generating accurate numerical simulations.

The authors present a detailed description of the elements necessary to analytically calculate biased power spectra with loop corrections. These analytical predictions are compared with the results of fully non-perturbative simulations, validating the approach. Finally, they propose a strategy to extend existing HEFT-based emulators for $\Lambda$CDM, such as \texttt{bacco} and \texttt{Aemulus}, to cosmologies beyond the standard model, opening the door to a more robust analysis of data from future galaxy surveys in the context of modified gravity theories.