We study level statistics in ensembles of integrable $N\times N$ matrices
linear in a real parameter $x$. The matrix $H(x)$ is considered integrable if
it has a prescribed number $n>1$ of linearly independent commuting partners
$H^i(x)$ (integrals of motion) $\left[H(x),H^i(x)\right] = 0$, $\left[H^i(x),
H^j(x)\right]$ = 0, for all $x$. In a recent work, we developed a
basis-independent construction of $H(x)$ for any $n$ from which we derived the
probability density function, thereby determining how to choose a typical
integrable matrix from the ensemble. Here, we find that typical integrable
matrices have Poisson statistics in the $N\to\infty$ limit provided $n$ scales
at least as $\log{N}$; otherwise, they exhibit level repulsion. Exceptions to
the Poisson case occur at isolated coupling values $x=x_0$ or when correlations
are introduced between typically independent matrix parameters. However, level
statistics cross over to Poisson at $ \mathcal{O}(N^{-0.5})$ deviations from
these exceptions, indicating that non-Poissonian statistics characterize only
subsets of measure zero in the parameter space. Furthermore, we present strong
numerical evidence that ensembles of integrable matrices are stationary and
ergodic with respect to nearest neighbor level statistics.Comment: 18 pages, 26 figures, discussion on number variance added; published
versio