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 ₁,…,O _k on , and we average their action over the N-particles. We show that, for every fixed , expectations of products of functions of the averaged observables approach, as N→∞, 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 ₁,…,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.