  We consider the time evolution of $N$ bosons in the mean field regime for
factorized initial data. In the limit of large $N$, the many body evolution can
be approximated by the non-linear Hartree equation. In this paper we are
interested in the fluctuations around the Hartree dynamics. We choose $k$
self-adjoint one-particle operators $O_1, \dots, O_k$ on $L^2 (\R^3)$, and we
average their action over the $N$-particles. We show that, for every fixed $t
\in \R$, expectations of products of functions of the averaged observables
approach, as $N \to \infty$, expectations with respect to a complex Gaussian
measure, whose covariance matrix can be expressed in terms of a Bogoliubov
transformation describing the dynamics of quantum fluctuations around the mean
field Hartree evolution. If the operators $O_1, \dots, O_k$ commute, the
Gaussian measure is real and positive, and we recover a "classical"
multivariate central limit theorem. All our results give explicit bounds on the
rate of the convergence (we obtain therefore Berry-Ess{\'e}en type central
limit theorems).
