Estimation of the remainder of a cubature formula on a Chebyshev grid for two-variable functions

For a cubature formula of the form
$$ \int_0^{2\pi}\int_0^{2\pi}f(x,y)\,dx\,dy=\frac{4\pi^2}{mn}\sum_{i=0}^{n-1}\sum_{j=0}^{m-1} f\biggl(\frac{2\pi i}{n},\frac{2\pi j}{m}\biggr)+R_{n,m}(f). $$
on a Chebyshev grid, the remainder $R_{n,m}(f)$ is proved to satisfy the sharp estimate
$$ \sup_{f\in H(r_1,r_2)}|S_{n,m}(f)|=O(n^{-r_1+1}+m^{-r_1+1}) $$
in some class of functions $H(r_1,r_2)$ defined by a generalized shift operator. Here, $r_1,r_2>1$; $\lambda^{-1}\le n/m\le\lambda$ with $\lambda>0$ and the constant in the $O$-term depends only on $\lambda$.