  We study global fluctuations for singular values of $M$-fold products of
several right-unitarily invariant $N \times N$ random matrix ensembles. As $N
\to \infty$, we show the fluctuations of their height functions converge to an
explicit Gaussian field, which is log-correlated for $M$ fixed and has a white
noise component for $M \to \infty$ jointly with $N$. Our technique centers on
the study of the multivariate Bessel generating functions of these spectral
measures, for which we prove a central limit theorem for global fluctuations
via certain conditions on the generating functions. We apply our approach to a
number of ensembles, including square roots of Wishart, Jacobi, and unitarily
invariant positive definite matrices with fixed spectrum, using a detailed
asymptotic analysis of multivariate Bessel functions to verify the necessary
conditions.
