Let us give a name to the partial sums
$$
S_{P,Q}(f)=\frac 1{PQ}\sum_{i=1}^P\sum_{j=1}^Q f(X_i,Y_j).
$$
and define the functions
$$
f_1\colon x\mapsto \mathbb E\left[f(x,Y_1)\right]-\mu, \quad, f_2\colon y\mapsto \mathbb E\left[f(X_1,y)\right]-\mu,
$$
$$
f_3(x,y)=f(x,y)-f_1(x)-f_2(y)-\mu
$$
(this look very similar to Hoeffding's decomposition for U-statistics).
We thus get
$$
S_{P,Q}(f)-\mu=\frac 1P\sum_{i=1}^Pf_1(X_i)+\frac 1Q\sum_{j=1}^Q f_2(Y_j)+S_{P,Q}(f_3).
$$
We can show that for each positive $\varepsilon$,
$$
\sum_{M,N\geqslant 1}\mathbb P\left(\frac 1{2^{M+N}}\max_{1\leq P\leq 2^M}\max_{1\leq Q\leq 2^N} \lvert S_{P,Q}(f_3)\rvert >\varepsilon\right)<\infty.
$$
Indeed, we first use Chebyshev's inequality to reduce this to prove that
$$\tag{*}
\sum_{M,N\geqslant 1}\mathbb E\left(\frac 1{2^{2(M+N)}}\max_{1\leq P\leq 2^M}\max_{1\leq Q\leq 2^N} \lvert S_{P,Q}(f_3)\rvert^2 \right)<\infty.
$$
To do so, observe that $\left(\max_{1\leq P\leq 2^M} \lvert S_{P,Q}(f_3)\rvert^2\right)_{Q\geqslant 1}=(M_Q)_{Q\geqslant 1}$ is a submartingale for the filtration
$\left(\mathcal F_Q\right)$, where $\mathcal F_Q=\sigma(X_i,1\leq i\leq 2^M,
Y_j,1\leqslant j\leqslant Q)$: by Jensen's inequality and monotonicity of conditional expectation,
$$
\mathbb E\left[M_Q\mid\mathcal F_{Q-1}\right]
\geqslant \max_{1\leq P\leq 2^M} \lvert \mathbb E\left[S_{P,Q}(f_3)\mid\mathcal F_{Q-1}\right]\rvert^2,
$$
and $\mathbb E\left[S_{P,Q}(f_3)\mid\mathcal F_{Q-1}\right]=S_{P,Q-1}(f_3)
+\sum_{i=1}^P \mathbb E\left[f_3(X_i,Y_{Q})\mid\mathcal F_{Q-1}\right]=S_{P,Q-1}(f_3)$, since $\mathbb E\left[Z\mid\mathcal F\vee\mathcal G\right]=\mathbb E\left[Z\mid\mathcal F\right]$ if $\mathcal G$ is independent of $\sigma(Z)\vee\mathcal F$,
$$
\mathbb E\left[f_3(X_i,Y_{Q})\mid\mathcal F_{Q-1}\right]=
\mathbb E\left[f_3(X_i,Y_{Q})\mid X_i\right]=0.
$$
Similarly, we can show that $\left(\lvert S_{P,2^N}(f_3)\rvert^2\right)_{P\geqslant 1}$ is a sub-martingale. By Doob's inequality, (*) reduces to
$$\tag{**}
\sum_{M,N\geqslant 1}\mathbb E\left(\frac 1{2^{2(M+N)}} \lvert S_{2^M,2^N}(f_3)\rvert^2 \right)<\infty.
$$
This follows from the fact that the summands are pairwise orthogonal.
In conclusion, the wanted convergence holds without imposing any rate on the convergence of $P$ and $Q$ to infinity. Some further remarks:
- The convergence for $f_3$ holds actually if only one of the indexes $P$ and $Q$ goes to infinity.
- Using Marcinkiewicz strong law of large numbers, we could show that for $1\leq p<2$,
$$
\min\{P^{1-1/p},Q^{1-1/p}\}\left(S_{P,Q}(f)-\mu\right)\to 0.
$$
- The convergence $S_{P,Q}(f)-\mu$ can hold under less restrictive assumptions than a moment of order two.