The notion of a $U$-statistic for an $n$-tuple of identical quantum systems is introduced in analogy to the classical (commutative) case: given a selfadjoint `kernel' $K$ acting on $(\mathbb{C}^{d})^{\otimes r}$ with $r<n$, we define the symmetric operator $U_{n}= {n \choose r} \sum_βK^{(β)}$ with $K^{(β)}$ being the kernel acting on the subset $β$ of $\{1,\dots ,n\}$. If the systems are prepared in the i.i.d state $ρ^{\otimes n}$ it is shown that the sequence of properly normalised $U$-statistics converges in moments to a linear combination of Hermite polynomials in canonical variables of a CCR algebra defined through the Quantum Central Limit Theorem. In the special cases of non-degenerate kernels and kernels of order $2$ it is shown that the convergence holds in the stronger distribution sense. Two types of applications in quantum statistics are described: testing beyond the two simple hypotheses scenario, and quantum metrology with interacting hamiltonians.