combined with sophisticated simulations in order to predict the spatial distribution of galaxies inside gravitationally collapsed dark matter halos [@0004-637X-575-2-587; @0004-637X-609-1-35; @0004-637X-633-2-791]. Furthermore, when MG configurations are considered, the fifth forces introduce an additional layer of complexity, scale-dependent growth is generally present even at the linear level and in the nonlinear scales, one needs to solve the scalar field Klein-Gordon (KG) equation that is highly nonlinear and adds to the computational costs significantly. The intermediate, quasi-linear, scales can fortunately be analytically accessed by higher order Perturbation Theory (PT) [@Bernardeau:2001qr; @2009PhRvD..80d3531C] approaches or hybrid methods [@Tassev:2013pn; @Valogiannis:2016ane].

The Lagrangian Perturbation Theory approach to structure formation is one the oldest and most popular analytical frameworks in the literature, that has been been particularly successful at describing the Baryon Acoustic Oscillation (BAO) peak [@Eisenstein:2005su], observed at a comoving scale of $\sim$110 Mpc/$h$, and the power-law correlation function, on comoving scales $\sim$20-90 Mpc/$h$, in $\Lambda$CDM [@Vlah:2014nta]. Combined with a model for halo bias [@1984ApJ...284L...9K; @Mo:1996cn; @Sheth:1999mn; @2009JCAP...08..020M; @PhysRevD.88.023515; @Matsubara:2008wx], it can be used to predict the 2-point statistics for halos in the real and redshift space [@2013MNRAS.429.1674C; @Matsubara:2008wx; @Vlah:2016bcl], which serves as a crucial step to the theoretical description of galaxy clustering. Additional contributions 