Sharp estimates for the convergence rate of “hyperbolic” partial sums of double fourier series in orthogonal polynomials

Two-variable functions $f(x,y)$ from the class $L_2=L_2((a,b)\times(c,d);p(x)q(y))$ with the weight $p(x)q(y)$ and the norm
$$ ||f||=\sqrt{\int_a^b\int_c^dp(x)q(x)f^2(x,y)dx\,dy} $$
are approximated by an orthonormal system of orthogonal $P_n(x)Q_n(y)$, $n, m=0, 1,\dots$, with weights $p(x)$ and $q(y)$. Let
$$ E_N(f)=\inf_{P_N}||f-P_N|| $$
denote the best approximation of $f\in L_2$ algebraic polynomials of the form
\begin{gather*} P_N(x,y)=\sum_{0<n,m<N}a_{m,n}x^ny^m,\\ P_1(x,y)=\mathrm{const}. \end{gather*}
Consider a double Fourier series of $f\in L_2$ in the polynomials $P_n(x)Q_m(y)$, $n, m=0, 1,\dots$, and its “hyperbolic” partial sums
\begin{gather*} S_1(f; x,y)=c_{0,0}(f)P_0(x)Q_0(y),\\ S_N(f; x,y)=\sum_{0<n,m<N}c_{n,m}(f)P_n(x)Q_m(y),\qquad N=2,3,\dots. \end{gather*}
A generalized shift operator $F_h$ and a $k$th-order generalized modulus of continuity $\Omega_k(A,h)$ of a function $f\in L_2$ are used to prove the following sharp estimate for the convergence rate of the approximation:
\begin{gather*} E_N(f)\leqslant(1-(1-h)^{2\sqrt{N}})^{-k}\,\Omega_k(f; h), \qquad h\in(0,1),\\ N=4,5,\dots;\qquad k=1,2,\dots. \end{gather*}
Moreover, for every fixed $N=4,9,16,\dots$, the constant on the right-hand side of this inequality is cannot be reduced.