Propagation of high-frequency wave in periodic media is a challenging problem due to the existence of multiscale characterized by short wavelength, small lattice constant and large physical domain size. Conventional computational methods lead to extremely expensive costs, especially in high dimensions. In this paper, based on Bloch decomposition and asymptotic analysis in the phase space, we derive the frozen Gaussian approximation for high-frequency wave propagation in periodic media and establish its converge to the true solution. The formulation leads to efficient numerical algorithms, which are presented in a companion paper [5].1 2 RICARDO DELGADILLO, JIANFENG LU, AND XU YANG avoid the boundary effects. Therefore the total number of grid points is huge, which usually leads to unaffordable computational cost, especially in high (d > 1) dimensions.An alternative efficient approach is to solve (1.1) asymptotically by the Bloch decomposition and modified WKB methods [3,4,7], which lead to eikonal and transport equations in the semi-classical regime. An advantage of this method is that the computational cost is independent of ε. However, the eikonal equation can develop singularities which make the method break down at caustics. The Gaussian beam method (GBM) [31] was then introduced by Popov to overcome this drawback at caustics. The idea is to allow the phase function to be complex and choose the imaginary part properly so that the solution has a Gaussian profile; see [15-20, 28, 37, 38] for recent developments. Similar ideas can be also found in the Hagedorn wave packet method [8,9]. Unlike the geometric optics based method, the Gaussian beam method allows for accurate computation of wave function around caustics [6,33]. But the problem is that the constructed beam must stay near the geometric rays to maintain accuracy. This becomes a drawback when the solution spreads [23,28,32].The Herman-Kluk propagator [11,21,22] was proposed for Schrödinger equation without the oscillatory periodic background potential. The method was rigorously analyzed in [35,36] and further extended as the frozen Gaussian approximation (FGA) for general high frequency wave propagation in [23][24][25]. The FGA method uses Gaussian functions with fixed widths, instead of using those that might spread over time, to approximate the wave solution. Despite its superficial similarity with the Gaussian beam method, it is different at a fundamental level. FGA is based on phase plane analysis, while GBM is based on the asymptotic solution to a wave equation with Gaussian initial data. In FGA, the solution to the wave equation is approximated by a superposition of Gaussian functions living in phase space, and each function is not necessarily an asymptotic solution, while GBM uses Gaussian functions (called beams) in physical space, with each individual beam being an asymptotic solution to the wave equation. The main advantage of FGA over GBM is that the problem of beam spreading no longer exists.In this paper, we extend FGA for computation of hi...