  Recent experimental progress using ultracold gases in optical lattices
necessitates a quantitative theoretical description for a significant number of
bosons. In the present paper, we investigate if time-dependent semiclassical
initial value methods, with propagators expressed as integrals over phase space
using classical trajectories, is suitable to describe interacting bosons,
concentrating on a single mode. Despite the nonlinear contribution from the
self-interaction, the corresponding classical dynamics allows for a largely
analytical treatment of the semiclassical propagator. We find that the
Herman-Kluk (HK) propagator conserves unitarity in the semiclassical limit
($n\to\infty$), but a decay of the norm is seen for low $n$. The frozen
Gaussian approximation (FGA, i.e. HK with unit prefactor) is explicitly shown
to violate unitarity in the present system for non-vanishing interaction
strength, even in the semiclassical limit. Furthermore, we show by evaluating
the phase space integral in steepest descent approximation, that the HK
propagator reproduces the exact spectrum correctly in the limit $n\to \infty$.
An error is, however, incurred in next-to-next-to-leading order (small
parameter $1/n$), as seen upon numerical evaluation of the integral or by
considering finite $n$ corrections to steepest descent. The FGA, in contrast,
is only accurate to lowest order, and an erroneous next-to-leading order term
in the energy spectrum was found analytically. Finally, as an example
application, we study the dynamics of wave packets by computing the time
evolution of the Wigner function. While the often-used truncated Wigner
approximation cannot capture any interferences present in the exact quantum
mechanical solution (known analytically), we find that the HK approach, despite
also using classical information only, reproduces the salient features of the
exact solution correctly.
