In this article we obtain a superexponential rate of convergence in total variation between the traces of the first m powers of a n x n random unitary matrices and a 2m-dimensional Gaussian random variable. This generalizes previous results in the scalar case to the multivariate setting, and we also give the precise dependence on the dimensions m and n in the estimates with explicit constants. We are especially interested in the regime where m grows with n and our main result basically states that if m <<root n, then the rate of convergence in the Gaussian approximation is Gamma(n/m + 1)(-1) times a correction. We also show that the Gaussian approximation remains valid for all m << n(2/3) without a fast rate of convergence.