  This work is concerned with asymptotic approximations of the semi-classical
Schr\"odinger equation in periodic media using Gaussian beams. For the
underlying equation, subject to a highly oscillatory initial data, a hybrid of
the Gaussian beam approximation and homogenization leads to the Bloch
eigenvalue problem and associated evolution equations for Gaussian beam
components in each Bloch band. We formulate a superposition of Bloch-band based
Gaussian beams to generate high frequency approximate solutions to the original
wave field. For initial data of a sum of finite number of band eigen-functions,
we prove that the first-order Gaussian beam superposition converges to the
original wave field at a rate of $\epsilon^{1/2}$, with $\epsilon$ the
semiclassically scaled constant, as long as the initial data for Gaussian beam
components in each band are prepared with same order of error or smaller. For a
natural choice of initial approximation, a rate of $\epsilon^{1/2}$ of initial
error is verified.
