  We study the statistics of matrix elements of local operators in the basis of
energy eigenstates in a paradigmatic integrable many-particle quantum theory,
the Lieb-Liniger model of bosons with repulsive delta-function interaction.
Using methods of quantum integrability we determine the scaling of matrix
elements with system size. As a consequence of the extensive number of
conservation laws the structure of matrix elements is fundamentally different
from, and much more intricate than, the predictions of the eigenstate
thermalization hypothesis for generic models. We uncover an interesting
connection between this structure for local operators in interacting integrable
models, and the one for local operators that are not local with respect to the
elementary excitations in free theories. We find that typical off-diagonal
matrix elements $\langle\boldsymbol{\mu}|O|\boldsymbol{\lambda}\rangle$ in the
same macro-state scale as $\exp(-c^{
O}L\ln(L)-LM^{O}_{\boldsymbol{\mu},\boldsymbol{\lambda}})$ where the
probability distribution function for
$M^{O}_{\boldsymbol{\mu},\boldsymbol{\lambda}}$ are well described by Fr\'echet
distributions and $c^{O}$ depends only on macro-state information. In contrast,
typical off-diagonal matrix elements between two different macro-states scale
as $\exp(-d^{ O}L^2)$, where $d^{O}$ depends only on macro-state information.
Diagonal matrix elements depend only on macro-state information up to
finite-size corrections.
