Summary: We develop an essentially optimal numerical method for solving two-scale Maxwell wave equations in a domain \(D \subset \mathbb{R}^d\). The problems depend on two scales: one macroscopic scale and one microscopic scale. Solving the macroscopic two-scale homogenized problem, we obtain the desired macroscopic and microscopic information. This problem depends on two variables in \(\mathbb{R}^d\), one for each scale that the original two-scale equation depends on, and is thus posed in a high dimensional tensorized domain. The straightforward full tensor product finite element (FE) method is exceedingly expensive. We develop the sparse tensor product FEs that solve this two-scale homogenized problem with essentially optimal number of degrees of freedom, i.e. the number of degrees of freedom differs by only a logarithmic multiplying factor from that required for solving a macroscopic problem in a domain in \(\mathbb{R}^d\) only, for obtaining a required level of accuracy. Numerical correctors are constructed from the FE solution. We derive a rate of convergence for the numerical corrector in terms of the microscopic scale and the FE mesh width. Numerical examples confirm our analysis.