High order harmonic generation in semiconductors
driven at near- and mid-IR wavelengths

The SC-TDDFT method considers the interaction of a bulk crystal unit cell with a spatially uniform electric field Bertsch et al. (2000); Otobe et al. (2008). This method has been used in the study of nonlinear and ultrafast dynamics Floss et al. (2018, 2019); Otobe (2012, 2016); Tancogne-Dejean et al. (2017); Yamada and Yabana (2020); Uemoto et al. (2019); Yamada and Yabana (2019); Schultze et al. (2014); Wachter et al. (2014). Utilizing the dipole approximation, the electronic motion can be described through the Bloch orbitals of the bulk, $u_{n\mathbb{k}}(\mathbb{r},t)$, where the wave vector $\bm{k}$ is contained within the three dimensional (3D) Brillouin zone of the unit cell. The TDKS equation for such orbitals is written as:

The applied electric field in this approach is defined by the vector potential of the incident pulse, $A(t)$. The electric potential $\phi(\mathbb{r},t)$ includes the Hartree potential from the electrons and the local contribution of the ionic potential. The terms $\delta V_{\text{ion}}$ and $V_{xc}(r,t)$ describe the nonlocal part of the ionic pseudopotential Troullier and Martins (1991) and the exchange-correlation potential Perdew and Zunger (1981), respectively. In our calculations, this exchange-correlation potential utilizes the adiabatic local-density approximation (LDA) Onida et al. (2002). Note that the LDA underestimates the bandgap of Si, as it does for diamond and other dielectrics. The experimental direct bandgap of Si is 3.4 eV while LDA gives 2.4 eV. We neglect the exchange-correlation term in the vector potential $A_{xc}(t)$ for simplicity and do not rigorously model the exchange-correlation effects of these types for infinite periodic systems Ullrich (2012); Vignale and Kohn (1996).

The averaged electric current density, used in generating HHG spectra, is calculated as follows:

Here $\Omega$ is the volume of the unit cell and $n$ is the band index restricted to occupied bands. The term $\delta J(t)$ denotes the contribution to the current density from the nonlocal part of the pseudopotential Bertsch et al. (2000) and is given as follows:

We utilize the open-source software package SALMON (Scalable Ab initio Light-Matter simulator for Optics and Nanoscience) Noda et al. (2019). In this code, a uniform 3D spatial grid is defined for the electron orbitals. The time evolution of the electron orbitals is carried out using the Taylor expansion method Yabana and Bertsch (1996).

For TDDFT calculations without dephasing, we define a cubic unit cell of the diamond structure containing eight atoms with varying side lengths depending on the specific target. For Di and Si calculations, the unit cell dimensions are defined as $a_{\text{Di}}=0.357$ nm and $a_{\text{Si}}=0.543$ nm. A spatial grid of $16^{3}$ points is used in Si and $24^{3}$ in Di calculations to enable an accurate description of the Bloch orbitals. Convergence over spatial (r-point) grids is observed at these grid sizes. Similarly, the 3D Brillouin zone is sampled by a grid made up of $32^{3}$ $k$-points. Note that convergence over the Brillouin zone is not observed at this grid size for diamond although qualitatively this size appears reasonable where we restrict ourselves to below the $50^{\text{th}}$ order harmonic. Otobe (2012) finds similar grid parameters for diamond but with convergence observed over the Brillouin zone. The timestep is set for stability of the calculation to be $2.5\times 10^{-3}$ fs for Si and $1\times 10^{-3}$ fs for Di. Note that tests for convergence are important in these calculations due to the discretization of the spatial, $k$-space and time parameters discussed here.

For the Si calculation with dephasing, a supercell of size $4^{3}$ describing a $512$ atom system is used. The Brillouin zone is sampled by a grid made up of $8^{3}$ $k$-points. Atomic positions are shifted from the equilibrium position with the diamond structure while are frozen during the time evolution of electron orbitals. They are generated from molecular dynamics simulations. The detail of the preparation will be described below.

In our SC-TDDFT calculations, we utilize the following incident pulse profile:

Here $E_{0}$ is the peak amplitude of the electric field, $\omega$ is the carrier frequency and $T$ is the full duration of the pulse. For the $\cos^{6}$ envelope, the FWHM duration $T_{\rm FWHM}$ is related to the full duration T by $T_{\rm FWHM}\simeq 0.3T$. In our calculations, we consider the photon energies $\hbar\omega=0.43$ eV ($\lambda=3000$ nm), $\hbar\omega=0.62$ eV ($\lambda=2000$ nm) or $\hbar\omega=1.55$ eV ($\lambda=800$ nm). We select typical pulse durations $T\gtrsim 100$ fs which are considerably larger than in previous HHG studies ($\sim 32$ fs in Floss et al. (2018, 2019) and 30 fs in Yamada and Yabana (2021)).

In the SC-TDDFT formalism, the HHG spectrum is calculated as a Fourier transform of the matter current

where $w(t)$ is the window smoothing function applied to remove spurious peaks in the spectra. A Gaussian convolution of the form

is applied to carry out frequency averaging as often the case in experimental HHG measurements.

We first consider bulk Si driven by a 2000 nm laser pulse with 200 fs full duration with a peak intensity of $\rm{5}\times 10^{{11}}~{}W/cm^{2}$. In  Fig. 1, the time profile of the vector potential is shown in the upper panel, and the induced matter current is shown in the lower panel. Here the simulated matter current displays large noise contributions that become particularly prominent after around $110$ fs. Note also that the current is very asymmetric compared to the symmetric driving pulse centered around $100$ fs. We thus find that the SC-TDDFT method assuming a diamond structure without any distortion effects records unphysically oscillating current towards the end of the specified pulse length. In Sec. III.C, we will show that the oscillating current mostly disappears once we introduce the dephasing effect caused by the thermal motion of atoms.

Fig. 2 displays the corresponding HHG spectrum (in green) generated from the matter current shown in Fig. 1 using Eqs. (5-7). We confirmed the simulation’s convergence over spatial (r-point) and momentum (k-point) grids up to $51^{\text{st}}$ order harmonic for Si at 2000 nm. Nevertheless we show the spectrum up to $99^{\text{th}}$ order to demonstrate both the presence of, and our methods application to, higher order harmonics. Although HHG signals can be identified in the unaltered spectrum (green) even around the highest order shown in the figure, the harmonic structure appears noisy and some harmonics cannot be clearly seen in this spectrum. Gaussian convolution using Eq. (8) overcomes this problem. Figure 2 shows the spectrum with the Gaussian filter (in red).

Fig. 3 demonstrates the major effect that pulse duration has on the clarity and extent of HHG spectra, considering a $3000$ nm pulse interacting with Si. The $200$ fs calculation (in green) shows clear harmonic peaks beyond the $27^{\text{th}}$ harmonic order where the $70$ fs calculation generally shows unclear harmonics until the $15^{\text{th}}$ order. In addition to this, we note that the harmonic cutoff, comparing the simulations for the $200$ fs and $70$ fs pulses, appears dependent on pulse duration. We further examined this dependence through the application of an inverse Gabor transform focused to reveal the underlying current contributions associated with the highest order harmonics. Our analysis demonstrates, as expected, the importance of multi-cycle current accumulation in the formation of harmonic structure, where we observe that a $70$ fs, $3000$ nm pulse possesses too few cycles to capture the higher order returns.

Our Si results at 2000 nm (Fig. 2) and 3000 nm (Fig. 3) need to be contrasted with the earlier calculation on Si at 3000 nm presented by Tancogne-Dejean et al. (2017). In particular, using equivalent parameters to those reported in Tancogne-Dejean et al. (2017), a direct comparison is achieved for 3000 nm in Fig. 3. Tancogne-Dejean et al. (2017) recorded noisy harmonics of the lower order while their higher order harmonics appeared clean. They related this effect with the joint density of states (JDOS) as a measure of the overlap between the valence and conduction bands. A large overlap and high JDOS would promote the recombination and inter-band HHG. Somewhat counter-intuitively, Tancogne-Dejean et al. (2017) discovered that a high JDOS was in fact detrimental to HHG emission making the corresponding harmonic peaks noisy. Conversely, those harmonics falling into a low JDOS region of the HHG spectrum stayed clean. This observation led the authors to conclude that it was the intra-band mechanism that was largely responsible for HHG in silicon driven at MIR. However, our HHG spectra at both $2000$ nm (Fig. 2) and $3000$ nm (Fig. 3) at $200$ fs total duration are both entirely clean even without averaging. Note that we do not observe a clear and obvious JDOS effect in our calculations at both 2000 nm and 3000 nm when the pulse duration is sufficiently large.

In Fig. 4 we display the HHG spectra of Di produced by the SC-TDDFT. The primary calculation is driven by an $800$ nm incident pulse of a $200$ fs full duration and the peak intensity of $\rm{1}\times 10^{{13}}~{}W/cm^{2}$. The Gaussian convoluted spectrum again reveals clear harmonics up to $41^{\text{st}}$ order of the fundamental frequency $\omega$, shown in the top panel of Fig. 4. The lower order harmonics of this spectrum are zoomed in on the bottom panel. The HHG spectrum is compared with an analogous calculation of Floss et al. (2019). The three sets of the SC-TDDFT calculations displayed in Fig. 4 are clearly differentiated by the pulse duration.

These authors employed the SC-TDDFT driven by a $\rm{2}\times 10^{{13}}~{}W/cm^{2}$ pulse with a total duration of $32$ fs Floss et al. (2018). Such a calculation did not produce a particularly clear HHG spectrum as is seen in Fig. 4. It is only the Maxwell propagation technique coupled with the SBE that allowed Floss et al. (2018) to generate a clear spectrum. We thus conclude that, with refined parameters and a longer pulse duration, one can obtain distinguishable high harmonics within the SC-TDDFT alone, while in previous calculations Floss et al. (2019) the Maxwell propagation or an inclusion of the explicit dephasing effect was required to obtain clear spectra that are often observed in measurements.

We expand the SC-TDDFT method to consider the atomic system at a finite temperature by allowing the atomic positions to be displaced from equilibrium. For this purpose, we consider bulk Si and first carry out a molecular dynamics simulation with an empirical force field at a specified temperature using a thermostat. We use LAMMPS software Plimpton (1995) to carry out the molecular dynamics simulation where a three-body Stillinger-Weber potential Stillinger and Weber (1985) is used. We pick up several atomic configurations with sufficiently long intervals, and use these configurations in the SC-TDDFT. This represents a first-principles approach introducing the dephasing effect into TDDFT for the first time without relying on phenomenological parameter $T_{2}$. This is a convenient approach where we consider that previously ultrafast dephasing times $T_{2}\sim 1$ fs had to be introduced to reconcile theoretical treatments and experiment. During the time evolution stage in the SC-TDDFT, the atomic positions are fixed. We carry out calculations using supercells of different sizes, and have found that results of different configurations are quite similar to each other if we choose sufficiently large supercell. We will show below the calculation results using a supercell containing 512 atoms. Employing $8^{3}$ $k$-points, the total number of atoms is equal to $262,144$, that is equal to the total number of atoms in a cell containing 8 atoms with the $32^{3}$ $k$-points that we showed in Fig. 2. In fact, we numerically confirmed that calculations without any distortion coincide with each other accurately between the two unit cells, 8 atoms with $32^{3}$ $k$-points and 512 atoms with $8^{3}$ $k$-points.

We show electron current density in bulk Si under the applied electric field in Fig. 5, caused by the same applied pulse as shown in Fig. 1. The black (blue) calculation displays the current under the atomic configuration at a temperature of 300K (180K). The matter current in disordered Si with distortion exhibits the desired dephasing after the pulse has ceased interacting. With dephasing present, signals beyond $150$ fs appear significantly suppressed and the noise contributions seen in Fig. 1 (red in Fig. 5) are removed.

Figure 6 demonstrates the effect of incorporating dephasing into the HHG calculations in Si where Gaussian convolutions are also applied. The black calculation shows the HHG spectrum using atomic positions from a molecular dynamics calculation at 300K. While signals up to $100^{\text{th}}$ order can be seen in Fig. 2 without atomic distortion, here we can see HHG signals up to $31^{\text{st}}$ order, and higher order signals cannot be identified. We also note that the HHG signals themselves become deeper between peaks and cleaner than those without atomic distortion, even when not incorporating a Gaussian convolution. Therefore, the dephasing effect modified the HHG spectrum in two ways, removing signals of HHG signals higher than $31^{\text{st}}$ order, and increasing peak-to-peak depth for those signals that survive.

In Fig. 6, the blue calculation shows the HHG spectrum using atomic positions from a molecular dynamics simulation at 180 K. Now we find HHG signals up to the $49^{\text{th}}$ order, where higher order signals are suppressed. This result at $180$ K clearly indicates that it is important to keep the material cold in measuring HHG to observe higher order signals. We however note that the distortion of atomic positions is generated by the classical molecular dynamics calculation. At very low temperatures, quantum zero-point motion will contribute significantly and should be included in discussing the HHG spectra.

The results with the effective thermally-induced dephasing effect, incorporated through atomic positions generated by molecular dynamics simulations, suggest that sample cooling could be used to greatly improve the high-order harmonic signal intensity observed in experiment. In principle, this should not pose a problem to current experiment as liquid nitrogen or helium cooling is commonly utilized in a variety of condensed matter physics contexts.

It should be noted that we have opted to consider fixed ionic coordinates in our dephasing simulations instead of performing evolution through classical mechanics. We do not incorporate an Ehrenfest molecular dynamics approach into TDDFT as we consider that spurious coherence may occur in the classical approximation when considering solids as, when one considers molecules, the coherence between two states will soon disappear as the ionic wave packets evolve differently between ionized and non-ionized electronic states and the spurious coherence that may occur in molecular simulations could extend to simulations concerning solids.