  We present a precise computation of the topological susceptibility
$\chi_{_\mathrm{YM}}$ of SU$(N)$ Yang-Mills theory in the large $N$ limit. The
computation is done on the lattice, using high-statistics Monte Carlo
simulations with $N=3, 4, 5, 6$ and three different lattice spacings. Two major
improvements make it possible to go to finer lattice spacing and larger $N$
compared to previous works. First, the topological charge is implemented
through the gradient flow definition; and second, open boundary conditions in
the time direction are employed in order to avoid the freezing of the
topological charge. The results allow us to extrapolate the dimensionless
quantity $t_0^2\chi_{_\mathrm{YM}}$ to the continuum and large $N$ limits with
confidence. The accuracy of the final result represents a new quality in the
verification of large $N$ scaling.
