= Solution
For a probability <measure-preserving system>, \b[<weak mixing> is the vanishing of averaged absolute correlation discrepancies]:
$$
\boxed{\frac1N\sum_{n=0}^{N-1}
\left|\mu(T^{-n}A\cap B)-\mu(A)\mu(B)\right|
\longrightarrow0
\quad(A,B\in\mathcal B).}
$$
By approximation with simple functions and the <Cauchy-Schwarz inequality>, this is equivalent to
$$
\frac1N\sum_{n=0}^{N-1}
\left|\langle U^nf,g\rangle-\left(\int f\,d\mu\right)
\overline{\left(\int g\,d\mu\right)}\right|\longrightarrow0
\qquad(f,g\in L^2).
$$
Here $\langle f,g\rangle=\int f\overline g\,d\mu$, and $U=U_T$. Absolute values are part of the definition: signed <Cesaro convergence of a sequence> of these correlations alone expresses <ergodicity>, not <weak mixing>.
Suppose the product system is <ergodic>. Fix a mean-zero $f\in L^2(\mu)$ and any $g\in L^2(\mu)$. On the product take
$$
F(x,y)=f(x)\overline{f(y)},\qquad
G(x,y)=g(x)\overline{g(y)}.
$$
These belong to $L^2(\mu\otimes\mu)$, and $\int F=|\int f|^2=0$. The <mean ergodic theorem> on the product, followed by pairing with $G$, gives
$$
\frac1N\sum_{n=0}^{N-1}|\langle U^nf,g\rangle|^2
=\left\langle\frac1N\sum_{n=0}^{N-1}(U\otimes U)^nF,G\right\rangle
\longrightarrow0.
$$
The <Cauchy-Schwarz inequality> in $n$ bounds the averaged absolute correlation by the square root of this quantity. Subtract the mean from a general $f$ to obtain the definition above. This proves \b[product <ergodicity> implies <weak mixing>], through the <square-correlation proof of weak mixing from product ergodicity>.
For the multiple averages, write $m(f)=\int f\,d\mu$ and use the following <Van der Corput lemma (Hilbert space sequences)>. If $(v_n)$ is bounded in a <Hilbert space>, every correlation average
$$
L_h=\lim_{N\to\infty}\frac1N\sum_{n=0}^{N-1}
\langle v_{n+h},v_n\rangle
$$
exists, and $H^{-1}\sum_{h=1}^H|L_h|\to0$, then $N^{-1}\sum_{n<N}v_n\to0$ in norm. One way to see the estimate is to replace $v_n$ by $H^{-1}\sum_{j=0}^{H-1}v_{n+j}$; for fixed $H$ the change in its long average tends to zero. The <Cauchy-Schwarz inequality> and expansion of the squared block norm give
$$
\limsup_{N\to\infty}\left\|\frac1N\sum_{n<N}v_n\right\|^2
\leq\frac{M^2}{H}
+\frac2H\sum_{h=1}^{H-1}\left(1-\frac hH\right)|L_h|,
\qquad M=\sup_n\|v_n\|.
$$
Let $H\to\infty$. This proves the auxiliary implication needed here.
A <weak mixing> system is <ergodic>: a mean-zero invariant $f$ would have the nonvanishing correlation $\langle U^nf,f\rangle=\|f\|_2^2$. Every positive power $T^r$ is also <weak mixing>, since for nonnegative correlation discrepancies $b_n$,
$$
\frac1H\sum_{h=0}^{H-1}b_{rh}
\leq\frac r{rH}\sum_{n=0}^{rH-1}b_n\longrightarrow0.
$$
This part of the <stability of weak mixing under powers and products> will control the induction.
In fact the stronger <arithmetic-progression multiple averages under weak mixing> holds:
$$
\boxed{\frac1N\sum_{n=0}^{N-1}
\prod_{j=1}^k U^{jn}f_j
\longrightarrow\prod_{j=1}^k m(f_j)
\quad\text{in }L^2,\qquad f_j\in L^\infty.}
$$
For $k=1$ this is the <mean ergodic theorem> and <ergodicity>. Suppose it is known for $k-1$ and first assume $m(f_k)=0$. Put $v_n=\prod_{j=1}^kU^{jn}f_j$. For fixed $h$, set
$$
g_{j,h}=(U^{jh}f_j)\overline{f_j}.
$$
Using invariance of the integral to remove the common $U^n$ gives
$$
\langle v_{n+h},v_n\rangle
=\int g_{1,h}\prod_{j=2}^k U^{(j-1)n}g_{j,h}\,d\mu.
$$
This identity does not require an inverse of $T$. Apply the induction hypothesis to the $k-1$ factors and pair their $L^2$ limit with $g_{1,h}$. It follows that
$$
L_h=\prod_{j=1}^k m(g_{j,h}),\qquad
|L_h|\leq
\left(\prod_{j=1}^{k-1}\|f_j\|_\infty^2\right)
|\langle U^{kh}f_k,f_k\rangle|.
$$
The average in $h$ of the right side tends to zero because $T^k$ is <weak mixing> and $m(f_k)=0$. The <Van der Corput lemma (Hilbert space sequences)> proves the zero $L^2$ limit. For general $f_k$, split it into $f_k-m(f_k)$ and its constant mean; the first term has the zero limit just proved, and the other term is $m(f_k)$ times the induction average for $k-1$. This completes the induction.
Pair this result with the real bounded $f_0$. For
$$
a_n=\int f_0\prod_{j=1}^kU^{jn}f_j\,d\mu,\qquad
\ell=\prod_{j=0}^km(f_j),
$$
\b[the requested <Cesaro limit> is]
$$
\boxed{\operatorname{C-lim}_{n\to\infty}a_n=\ell.}
$$
To obtain <convergence in density of a sequence>, we also need the product system to be <weak mixing>. For simple tensors, its correlations are products of single-system correlations. If $c_n\to c$ and $d_n\to d$ in averaged absolute discrepancy and both sequences are bounded, then
$$
|c_nd_n-cd|\leq |c_n-c|\,|d_n|+|c|\,|d_n-d|,
$$
so the product discrepancy has zero average. Finite sums of tensors are dense in $L^2(\mu\otimes\mu)$, and the <Cauchy-Schwarz inequality> extends the conclusion to arbitrary $L^2$ functions. Hence $T\times T$ is <weak mixing>.
Apply the multiple-average result to $F_j=f_j\otimes f_j$ on that product. Because the original $f_j$ are real,
$$
\frac1N\sum_{n<N}a_n^2
\longrightarrow\prod_{j=0}^k\left(\int F_j\,d(\mu\otimes\mu)\right)
=\ell^2.
$$
Together with the first-moment limit this gives
$$
\frac1N\sum_{n<N}|a_n-\ell|^2\longrightarrow0.
$$
The <mean-square criterion for convergence in density> now yields, for every $\varepsilon>0$,
$$
\frac1N\#\{0\leq n<N:|a_n-\ell|\geq\varepsilon\}
\leq\frac1{\varepsilon^2N}\sum_{n<N}|a_n-\ell|^2\longrightarrow0.
$$
Therefore \b[the <density convergence of multiple weak-mixing correlations> gives the same answer]:
$$
\boxed{\operatorname{D-lim}_{n\to\infty}a_n
=\prod_{j=0}^k\int f_j\,d\mu.}
$$
The second-moment argument is essential; the signed <Cesaro limit> alone would not imply this conclusion.
Back to article page