# Interference-Induced Phenomena in High-Order Harmonic Generation from Bulk Solids
## 2 Model and Method
### Numerical Approach  
In order to simplify the calculations, we are going to use a one-dimensional model in the following (see Figure 1). Note that this can be adequate for linearly polarized excitations. We consider two model potentials that have one point in common: both produce a band scheme with a (maximal) gap of \(3.2\,\mathrm{eV}\), which corresponds to the case of the ZnO crystal. The explicit form of the first potential is the following:  
\[U^{(a)}(x)=-U^{(a)}_{0}\sum_{i=1}^{2}\cos^{2}\left[\pi\left(x-x_{i}^{(a)} \right)/\Delta^{a}\right], \tag{6}\]  
where the functions \(\cos^{2}\) are assumed to be zero when their arguments are not in the interval \([-\pi/2,\pi/2]\). The parameters are: \(U^{(a)}_{0}=25\) eV, \(\Delta^{(a)}/a=0.15\), \(x_{1}^{(a)}/a=0.3\) and \(x_{2}^{(a)}/a=0.607\), where \(a\) denotes the lattice constant. \(U^{(b)}\) is also a double-well potential, as given by:  
\[U^{(b)}(x)=\mathcal{N}\left[-U^{(b)}_{shift}+U^{(b)}_{0}\sinh^{2}\left(\frac {x-x_{1}^{(b)}}{a}\right)\right]+\mathcal{N}\left[-U^{(b)}_{shift}+U^{(b)}_{0} \sinh^{2}\left(\frac{x-x_{2}^{(b)}}{a}\right)\right], \tag{7}\]where \(\mathcal{N}(x)\) denotes the negative part of \(x\), that is, \(\mathcal{N}(x)=x\) when \(x\leq 0\), and zero otherwise. Additionally, \(U_{shift}^{(b)}=93.861\) eV, \(U_{0}^{(b)}=4080.925\) eV, \(x_{1}^{(b)}/a=0.18\), and \(x_{2}^{(b)}/a=0.7\). The potentials, together with the band structure that they induced, are shown in Figure 2. Note that, apart from the bandgap, these potentials correspond to considerably different band schemes; thus, they can be used to check to what extent our results depend on the particular choice of the model. Inserting these potentials in the 1D version of Equation (2), we obtain the Bloch states \(\Psi_{n,k}(x)\). We use a finite number of \(k\) values, which are uniformly distributed in the first Brillouin zone \([-\pi/a,\pi/a]\), where \(a\) is the lattice constant. As a next step, we calculate the matrix elements \(\langle n,k|H_{0}|n^{\prime},k^{\prime}\rangle\) and  
Figure 1: Schematic view of an external laser-driven periodic structure. The duration of the pulse is denoted by \(\tau\) and the lattice constant of the 1D crystal is given by \(a\). Note that, usually, the central wavelength of the incoming laser pulse is larger than \(a\), meaning that dipole approximation can be used. Figure 2: Top panel: The model potentials \(U^{(a)}\) (red curve) and \(U^{(b)}\) (grey curve) as a function of the coordinate \(x\) in a unit cell. The corresponding band schemes are shown in the bottom panel (where an additive scaling factor has been used in order to see more clearly that the band gaps are the same). \(\langle n,k|p_{x}|n^{\prime},k^{\prime}\rangle\), which are both proportional to \(\delta_{kk^{\prime}}\). Using these matrix elements and the initial conditions (3), the dynamical equation (4) can be integrated by numerical means. In these calculations, we assume that the only nonzero component of the vector potential is given by:  
\[A_{x}(x,t)=\left\{ \begin{array}{l}A_{0}\sin^{2}(\pi t/\tau)\sin(\omega_{L}t+ \varphi),\,\mbox{if $t\in[0,\tau]$,}\\ 0\,\,\mbox{otherwise,}\end{array} \right., \tag{8}\]  
where \(\tau\) measures the duration of the pulse, and \(\omega_{L}=2\pi c/(3\ \upmu\mbox{m})\) denotes the central frequency of the excitation. In dipole approximation, the vector potential has no spatial dependence. Note that since we are considering relatively long pulses (\(\approx\)30 optical cycles at intensity FWHM), the results are expected to be almost independent of the carrier-envelope phase, \(\varphi\). Therefore, we use \(\varphi=0\) in the rest of the paper. When the time-dependent density matrix is obtained, the expectation value of the kinetic momentum can be calculated in a way that leads to the function \(J(t)\) (the 1D version of Equation (5)), and the fast Fourier transform provides the HHG spectra.