In this article we obtain a superexponential rate of convergence in total variation between the traces of the first m powers of a n × 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 (Formula Presented), then the rate of convergence in the Gaussian approximation is (Formula Presented) times a correction. We also show that the Gaussian approximation remains valid for all (Formula Presented) without a fast rate of convergence.
