for all quantities. The ensemble average K in (\[LX\]) can be cast into an exponent of a power series in cumulants, through the cumulant expansion theorem, $\langle e^{iX} \rangle = \exp\left[\sum_{N=1}^{\infty}\frac{i^N}{N!} \langle X^{N} \rangle_c \right]$, which, combined with a multinomial expansion, gives: $$\label{multinom}
K(\bold{q},\bold{k},\lambda_1,\lambda_2)=\exp \left[\sum \frac{i^{n+m+r}}{m!n!r!}\lambda_1^m\lambda_2^n k_{i_1}..k_{i_r}\langle \delta_1^{m} \delta_1^{n} \Delta_{i_1}..\Delta_{i_r}\rangle_c\right],$$ in terms of a series of correlators $$\begin{aligned}\label{eq:correl}
\sigma_R^2 &= \langle \delta^2 \rangle_c \\
\xi_L(\vec{q}) &= \langle \delta_1 \delta_2 \rangle_c,  \\
A_{ij}^{m n}(\vec{q}) &= \langle \delta_i^{m} \delta_j^{n} \Delta_i \Delta_j\rangle_c,  \\ 
W_{ijk}^{m n}(\vec{q}) &= \langle \delta_i^{m} \delta_j^{n} \Delta_i \Delta_j \Delta_k\rangle_c,  \\ 
U_{i}^{m n}(\vec{q}) &= \langle \delta_1^{m} \delta_2^{n} \Delta_i \rangle_c,
\end{aligned}$$ where we adopted the commonly used notation for the Lagrangian cumulants in (\[eq:correl\]).

Keeping all terms in (\[multinom\]) that contain cumulants up to third order, which is the equivalent of the one-loop correction to the linear power spectrum, results in a highly oscillatory integrand that presents challenges when ensuring the integral is fully converged. [@Matsubara:2007wj] proposed expanding all contributions to the exponent but the scale-independent “zero-lag" piece of $A_{ij}^{0 0}$, which results in a non-perturbative resummation scheme that is simpler to handle analytically. Building upon this result, [@2013MNRAS.429.1674C] proposed keeping all the terms of $A_{ij}^{0 0}$ in the exponent, in their 