In this article, we obtain a super-exponential rate of convergence in total variation between the traces of the first $m$ powers of an $n\times 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 estimate with explicit constants. We are especially interested in the regime where $m$ grows with $n$ and our main result basically states that if $m\ll \sqrt{n}$, then the rate of convergence in the Gaussian approximation is $Γ(\frac nm+1)^{-1}$ times a correction. We also show that the Gaussian approximation remains valid for all $m\ll n^{2/3}$ without a fast rate of convergence.