  We study the large-$N$ scaling behavior of the $\theta$ dependence of the
ground-state energy density $E(\theta)$ of four-dimensional (4D) $SU(N)$ gauge
theories and two-dimensional (2D) $CP^{N-1}$ models, where $\theta$ is the
parameter associated with the Lagrangian topological term. We consider its
$\theta$ expansion around $\theta=0$, $E(\theta)-E(0) = {1\over 2}\chi
\,\theta^2 ( 1 + b_2 \theta^2 + b_4\theta^4 +\cdots)$ where $\chi$ is the
topological susceptibility and $b_{2n}$ are dimensionless coefficients. We
focus on the first few coefficients $b_{2n}$, which parametrize the deviation
from a simple Gaussian distribution of the topological charge at $\theta=0$. We
present a numerical analysis of Monte Carlo simulations of 4D $SU(N)$ lattice
gauge theories for $N=3,\,4,\,6$ in the presence of an imaginary $\theta$ term.
The results provide a robust evidence of the large-$N$ behavior predicted by
standard large-$N$ scaling arguments, i.e. $b_{2n}= O(N^{-2n})$. In particular,
we obtain $b_2=\bar{b}_2/N^2 + O(1/N^4)$ with $\bar{b}_2=-0.23(3)$. We also
show that the large-$N$ scaling scenario applies to 2D $CP^{N-1}$ models as
well, by an analytic computation of the leading large-$N$ dependence.
