Hilbert series of this module, schematically of the form $$H(R_{\br}^{G_{\br}}/I_{\br},t) = \sum_{k=0}^{\infty}c_{k\, \br} t^k \sim \sum_{k=0}^{\infty}t^k \sum_{k_1+\dots+k_N = k}\big\{\text{conditions}\big\}.$$ Instead of resorting to involved mathematics, to obtain the coefficients $c_{k\, \br}$ we build on the experience gained from studying the previous simpler systems. Due to the EOM, the $\big(P_1^{(i)}\big)^{r_i+1}$ terms in the ideal require $k_i \le r_i$, as in eq. . As when we handled IBP alone, an independent IBP relation among the $\phi_1^{r_1} \cdots \phi_N^{r_N} \partial^k$ operators is generated from a total derivative acting on each of the operators with one less derivative. In this spirit, it is very tempting to write $$c_{k\, \textbf{r}}^{} \sim c_{k\, \textbf{r},\text{EOM}}^{} - c_{k-1\, \textbf{r},\text{EOM}}^{} =\sum_{k_1+\dots+k_N = k}I_{\{k_i\le r_i\}} - \sum_{k_1+\dots+k_N = k-1}I_{\{k_i\le r_i\}}.$$ However, it is almost immediately obvious that this expression cannot be correct, because it will go negative when $k$ is too large. When this happens, we should take $c_{k\, \textbf{r}}=0$. It turns out that the correction condition to guarantee $c_{k\, \br}$ be non-negative is $2k\le r\equiv r_1+\cdots +r_N$. Therefore, we have $$c_{k\, \textbf{r}}^{} \propto I_{2k \le r} \left( c_{k\, \textbf{r},\text{EOM}}^{} - c_{k-1\, \textbf{r},\text{EOM}}^{} \right).$$

Additionally, there is another constraint coming from the interplay of IBP and EOM. Using IBP, one 