# Interference-Induced Phenomena in High-Order Harmonic Generation from Bulk Solids
## 2 Model and Method
### The Source of the High Harmonic Signal  
In order to solve the dynamics induced by the Hamiltonian (1), the initial state has to specified. In other words, by expanding the quantum mechanical state \(\left|\Phi(t)\right\rangle=\Phi(\mathbf{r},t)\) of the electron in terms of the Bloch states \(\left|n,\mathbf{k}\right\rangle=\Psi_{n,\mathbf{k}}(\mathbf{r})\) as \(\left|\Phi(t)\right\rangle=\sum_{n,\mathbf{k}}c_{n,\mathbf{k}}(t)|n,\mathbf{k}\rangle\), the coefficients \(c_{n,\mathbf{k}}(t=0)\) should be given. (Note that, for the sake of simplicity, we assume periodic boundary conditions, and consequently the k-space will not be continuous--there will be a discrete, densely spaced series of k vectors.) However, the most plausible assumption for the initial state of the crystal is that it is in thermal equilibrium before the arrival of the laser pulse. Thermal states are not pure quantum mechanical states; thus, instead of a state vector, we have to use a density matrix, \(\rho\). Since the Bloch states are energy eigenstates of the crystal, all density matrices that represent thermal equilibria are diagonal on this basis. At room temperature, for a bandgap around 3 eV, the corresponding statistical weights are practically nonzero only for the valence bands (VB). That is, if we take a single VB into account with the corresponding band index \(n_{0}\), we can write:  
\[\rho(t=0)=\frac{1}{N}\sum_{\mathbf{k}}|n_{0},\mathbf{k}\rangle\langle n_{0}, \mathbf{k}|, \tag{3}\]  
where the constant \(N\) equals to the number of \(\mathbf{k}\) vectors in the first Brillouin zone, and provides normalization: \(\mathrm{Tr}\rho(0)=1\). Note that the term \"thermal equilibrium\" is not completely adequate for the state given by Equation (3). Strictly speaking, this state corresponds to the low temperature limit. However, for the large-bandgap case that we consider, the exact thermal state at room temperature is numerically the same as \(\rho(t=0)\) above. The time evolution of the density operator is governedby the von Neumann equation, supplemented by a phenomenological term that takes unavoidable decoherence effects into account:  
\[\frac{\partial}{\partial t}\rho(t)=-\frac{i}{\hbar}[H(t)\rho(t)-\rho(t)H(t)]+ \left.\frac{\partial}{\partial t}\rho(t)\right|_{dec}, \tag{4}\]  
where the second term on the right-hand side forces the population toward the valence band with a rate of \(\gamma_{d}\), and destroys the off-diagonal matrix elements (quantum mechanical coherences) with a rate of \(\gamma_{od}\). (Note that these effects do not change the index \(\mathbf{k}\) either.) For the calculations presented here, the diagonal and off-diagonal relaxation rates are \(\gamma_{d}=0.1\,\mathrm{fs}^{-1}\) and \(\gamma_{od}=0.3\,\mathrm{fs}^{-1}\), respectively. Note that \(\gamma_{od}\) is the inverse of the dephasing time \(T_{2}\), that can be estimated to fall within the range of 15-50 fs for low-energy electrons [47], but can be as short as a few fs, such as for the case of a strong longitudinal phonon-electron interaction in SiO\({}_{2}\)[48]. Although, according to our best knowledge, scattering rates for hot electrons have never been measured directly, \(T_{2}\approx 4\) fs have often been used in theoretical calculations, such as in [49; 50]. This is in agreement with our current choice. The diagonal relaxation described by \(\gamma_{d}\) describes energy relaxation, which is a slower effect that usually plays a less significant role on the fs timescale, but for the sake of completeness it is worth including it as well. As one can check, using the dipole approximation and velocity gauge, \(\rho\) will always be diagonal in the index \(\mathbf{k}\)--that is, the dynamics only mix the band indices (\"vertical transitions\"), but not different values of \(\mathbf{k}\). In a thermal equilibrium (without external bias), the sample does not radiate, and the net current flowing through it is zero.