  We employ stabilization methods and second order Poincar\'e inequalities to
establish rates of multivariate normal convergence for a large class of vectors
$(H_s^{(1)},...,H_s^{(m)})$, $s \geq 1$, of statistics of marked Poisson
processes on $\mathbb{R}^d$, $d \geq 2$, as the intensity parameter $s$ tends
to infinity. Our results are applicable whenever the constituent functionals
$H_s^{(i)}$, $i\in\{1,...,m\}$, are expressible as sums of exponentially
stabilizing score functions satisfying a moment condition. The rates are for
the $d_2$-, $d_3$-, and $d_{convex}$-distances. When we compare with a centered
Gaussian random vector, whose covariance matrix is given by the asymptotic
covariances, the rates are in general unimprovable and are governed by the rate
of convergence of $s^{-1} {\rm Cov}( H_s^{(i)}, H_s^{(j)})$,
$i,j\in\{1,...,m\}$, to the limiting covariance, shown to be of order
$s^{-1/d}$. We use the general results to deduce rates of multivariate normal
convergence for statistics arising in random graphs and topological data
analysis as well as for multivariate statistics used to test equality of
distributions. Some of our results hold for stabilizing functionals of Poisson
input on suitable metric spaces.
