  For smooth bounded domains in $\mathbb{R}$, we prove upper and lower $L^2$
bounds on the boundary data of Neumann eigenfunctions, and prove
quasi-orthogonality of this boundary data in a spectral window. The bounds are
tight in the sense that both are independent of eigenvalue; this is achieved by
working with an appropriate norm for boundary functions, which includes a
`spectral weight', that is, a function of the boundary Laplacian. This spectral
weight is chosen to cancel concentration at the boundary that can happen for
`whispering gallery' type eigenfunctions. These bounds are closely related to
wave equation estimates due to Tataru.
  Using this, we bound the distance from an arbitrary Helmholtz parameter $E>0$
to the nearest Neumann eigenvalue, in terms of boundary normal-derivative data
of a trial function $u$ solving the Helmholtz equation $(\Delta-E)u=0$. This
`inclusion bound' improves over previously known bounds by a factor of
$E^{5/6}$. It is analogous to a recently improved inclusion bound in the
Dirichlet case, due to the first two authors.
  Finally, we apply our theory to present an improved numerical implementation
of the method of particular solutions for computation of Neumann eigenpairs on
smooth planar domains. We show that the new inclusion bound improves the
relative accuracy in a computed Neumann eigenvalue (around the $42000$th) from
9 digits to 14 digits, with little extra effort.
