Ab-initio multi-scale simulation of high-harmonic generation in solids

TL;DR

This paper presents the first ab-initio multi-scale simulation of high-harmonic generation in solids, revealing how mesoscopic effects influence the formation of clean harmonic spectra in dense media.

Contribution

It introduces a novel multi-scale simulation approach for solid-state high-harmonic generation, emphasizing the importance of mesoscopic effects like Brillouin zone sampling and pulse propagation.

Findings

Mesoscopic effects significantly influence harmonic spectra.

Realistic Brillouin zone sampling is crucial.

Pulse propagation affects harmonic formation.

Abstract

High-harmonic generation by a highly non-linear interaction of infrared laser fields with matter allows for the generation of attosecond pulses in the XUV spectral regime. This process, well established for atoms, has been recently extended to the condensed phase. Remarkably well pronounced harmonics up to order ~30 have been observed for dielectrics. We present the first ab-initio multi-scale simulation of solid-state high-harmonic generation. We find that mesoscopic effects of the extended system, in particular the realistic sampling of the entire Brillouin zone, the pulse propagation in the dense medium, and the inhomogeneous illumination of the crystal have a strong effect on the formation of clean harmonic spectra. Our results provide a novel explanation for the formation of clean harmonics and have implications for a wide range of non-linear optical processes in dense media.

Equations

$${\partial }_t{\rho }_{mn}^{\vec{k}}$$

$$S_{\hat{n}}(\omega )\propto {\left|\mathfrak{Ft}\left\{\frac{d}{dt}\vec{J}(t)\cdot \hat{n}\right\}\right|}^2={\omega }^2{\left|\vec{J}(\omega )\cdot \hat{n}\right|}^2$$

$$\frac{1}{c^2}{\partial }_t^2\vec{A}(\vec{r},t)-{\nabla }^2\vec{A}(\vec{r},t)=\frac{4\pi }{c}\vec{J}(\vec{r},t)$$

$$\sigma (\omega )=\frac{\omega }{4\pi i}[\epsilon (\omega )-1].$$

$$J(t)=\frac{1}{2\pi }\int d\omega e^{-i\omega t}\sigma (\omega )F(\omega )$$

$$⟨\vec{J}(t)⟩=\frac{1}{R^2}{\int }_0^{\infty }d\rho \rho \vec{J}[t,I(\rho )]$$

Full text

The generation of high harmonics (HHG) in the non-linear interaction of intense ultrashort infrared (IR) laser pulses with matter has turned out to be a highly successful route towards the generation of attosecond pulses in the EUV and XUV spectral regimes Drescher et al. (2001); Zhao et al. (2012); Takahashi et al. (2013); Hammond et al. (2016). It has become the workhorse of investigation of a vast array of electronic processes on the attosecond time scale Krausz and Ivanov (2009). Expanding the range of accessible photon energies and intensities faces, however, fundamental limitations. Experimental and theoretical investigations have established a scaling of the cut-off energy $E_{\mathrm{cut}}\propto\lambda^{2}$ for HHG from atoms in the gas phase raising hopes to reach ever higher photon energies by increasing the wavelength $\lambda$ of the driving laser pulse. However, the intensity in the cut-off region was found to scale unfavorably $I_{\mathrm{cut}}\propto\lambda^{-5.3}$ due to the large spatial dispersion of the electron wave packet upon return to its parent atom Tate et al. (2007); Schiessl et al. (2007); Colosimo et al. (2008); Frolov et al. (2008); Ishikawa et al. (2009). Propagation effects in gas filled capillaries have been found to partially offset this suppression at high $\lambda$ Popmintchev et al. (2012).

Extending HHG to the condensed phase promises to overcome some of these limitations to enable compact and brighter light sources and to open up the novel field of solid-state photonics on the attosecond scale. The recent observation of HHG in solids for intensities below the damage threshold Ghimire et al. (2011); Mücke (2011); Ghimire et al. (2014); Schubert et al. (2014); Hohenleutner et al. (2015); Luu et al. (2015); Ndabashimiye et al. (2016) suggests opportunities for controlling electronic dynamics Hohenleutner et al. (2015); Luu et al. (2015) and for an all-optical reconstruction of the band structure Vampa et al. (2015a).

The observed solid-state HHG substantially differs from the corresponding atomic spectra. For example, while for atoms the cut-off frequency $\omega_{\mathrm{cut}}^{\mathrm{HHG}}$ scales linearly with the (peak) intensity $I_{0}$ of the driving pulse Krause et al. (1992); Corkum (1993), for HHG from bulk dielectrics or semiconductors $\omega_{\mathrm{cut}}^{\mathrm{HHG}}$ scales linearly with the peak field $F_{0}$ (or $\sqrt{I_{0}}$) Ghimire et al. (2014). The processes underlying solid-state HHG have remained a matter of debate. Several simplifying models have been proposed accounting for Bloch oscillations within a single band (“intraband harmonics”) Ignatov and Romanov (1976); Feise and Citrin (1999) and non-linear interband polarization (“interband harmonics”) Schubert et al. (2014); Golde et al. (2011); Vampa et al. (2015b) as sources of HHG. Most descriptions involve the semiconductor Bloch equations (SBE, Golde et al. (2008)) using input parameters on various levels of sophistication and a varying number of energy bands Hawkins et al. (2015); Wu et al. (2015). Recently, first simulations employing time-dependent density functional theory (TDDFT, Runge and Gross (1984)) have become available Otobe (2016); Tancogne-Dejean et al. (2017).

One major puzzle has remained so far unresolved: while many experiments display remarkably “clean” harmonic spectra with pronounced peaks near multiples of the driving frequency (odd multiples when inversion symmetry is preserved) all the way up to the cut-off frequency, corresponding simulations display a noisy spectrum lacking any clear harmonic structure over a wide range of frequencies in the “plateau” region above the band-gap energy. In previous theoretical works Vampa et al. (2014); Luu et al. (2015); Hohenleutner et al. (2015); Luu and Wörner (2016); Yu et al. (2016); You et al. (2017), this problem was addressed by proposing remarkably short dephasing times $T_{2}$ in the SBEs of the order of $T_{2}\approx 1$ fs or less than an optical half-cycle. While such short decoherence times yield “cleaner” harmonic spectra in qualitative agreement with the experiment, they raise important questions as to the ultrafast decoherence processes for electronic excitations in solids.

The point of departure of the present work is the first realistic ab-initio multi-scale simulation of HHG by self-consistently treating the microscopic non-linear response and the mesoscopic propagation of the optical signal. Diamond can serve as a prototypical bulk dielectric for which a full ab-initio treatment is feasible. On the microscopic scale we employ two complementary methods: an ab-initio TDDFT simulation and a multi-band SBE approach with input parameters from ab-initio ground-state DFT calculations and a dense 3-dimensional Brillouin zone (BZ) sampling. Properties of the extended solid target on the mesoscopic scale enter through propagation effects on the light fields in dense matter and the inhomogeneous field distribution within the focal spot that is typically smaller than the extended target. These mesoscopic effects turn out to have a surprisingly strong influence on the resulting HHG spectra.

We use TDDFT in a real-space real-time implementation Yabana et al. (2012); Wachter et al. (2015) to simulate the electronic dynamics driven by the strong IR field $\vec{F}(t)$ employing the adiabatic local-density approximation (LDA). Alternatively, we implement the SBEs by propagating the elements of the density matrix $\rho^{\vec{k}}_{mn}$,

[TABLE]

with the transition energy $\omega^{\vec{k}}_{mn}=\epsilon^{\vec{k}}_{m}-\epsilon^{\vec{k}}_{n}$ and the transition dipole elements $\vec{d}^{\;\vec{k}}_{mn}$ between Houston orbitals of the valence (VB) and conduction bands (CB), $n$, $|n,\vec{k}+\vec{A}/c\rangle$ Houston (1940); Krieger and Iafrate (1986) at wave vectors $\vec{k}$ displaced by the time-dependent vector potential $\vec{A}/c$. Both $\omega^{\vec{k}}_{mn}$ and $\vec{d}^{\;\vec{k}}_{mn}$ are taken from ground state (GS) ab-initio DFT serving also as initial state of the TDDFT calculation. The second term on the right-hand side describes the decoherence characterized by a dephasing time $T_{2}$ (a corresponding term accounting for the population relaxation with relaxation time $T_{1}$ is omitted because $T_{1}\gg T_{2}$).

The crystal is irradiated by pulses polarized along the $\Gamma-\mathrm{X}$ direction having a total duration of 32 fs (corresponding to $\tau_{\text{p}}\approx 6.8$ fs full-width at half maximum of the peak intensity), a carrier wavelength of 800 nm [photon energy $\omega_{\text{IR}}=1.55$ eV], peak intensities $I_{0}$ of up to $2\times 10^{13}$ W/cm2, and a $\sin^{6}$-envelope for the field. From the time-dependent induced current densities $\vec{J}(t)$ derived from either methods the resulting high harmonic spectrum is calculated using Larmor’s formula Jackson (1999)

[TABLE]

where $\hat{n}$ is the unit vector in polarization direction.

We account for the propagation through the extended system by combining the microscopic SBE solution for the current density $\vec{J}(t)$ of individual cells with the solution of Maxwell’s equations Yabana et al. (2012)

[TABLE]

where the cells are placed on a mesoscopic grid along the propagation direction [$\vec{A}(\vec{r},t)\rightarrow\vec{A}(X,t)$] with grid spacing $\Delta X=8$ nm and a crystal thickness of up to 1 $\mu$m. We use a 5-point stencil for the approximation of the second derivative in space and a standard 4-th order Runge-Kutta propagator to solve the differential equations in time. The microscopic responses to the impinging pulse at different grid points $i$, $\vec{J}(X_{i},t)$, are thus coupled via Eq. (3). In particular, high harmonic radiation emitted from the source term $\vec{J}(X_{i},t)$ propagates through the crystal via Eq. (3) driving the response at the other microscopic sites thereby accounting for re-absorption of HHG within the crystal. From the electric field $\vec{F}(\omega)$ of the transmitted pulse we retrieve the harmonic spectrum via $S_{\hat{n}}(\omega)\propto\omega^{2}|\vec{F}(\omega)\cdot\hat{n}|^{2}$.

Within a microscopic calculation employing a single cell of the periodic structure and neglecting propagation we have first verified that for moderate intensities of the driving pulse ($I_{0}\sim 10^{12}$ W/cm2) the SBE results for the time-dependent current density rapidly converge to that of the TDDFT prediction when consistent input for the band structure is used and dephasing is neglected ($T_{2}\rightarrow\infty$). Enforcing a continuous phase evolution of the dipole matrix elements along any trajectory through the BZ, in particular near narrow avoided crossings (see Fig. S1) and a fine $k$-grid ($\Delta k\simeq 0.01$) have been found to be crucial prerequisites for convergence to the ab-initio TDDFT simulation. Another important convergence parameter is the number of VBs and CBs Hawkins et al. (2015); Wu et al. (2015). The present results (Fig. 1) show, however, that approximations including only a one-dimensional cut through the BZ and, thus, without properly sampling the full BZ as have been frequently employed Hawkins et al. (2015); Luu et al. (2015); Luu and Wörner (2016); Yu et al. (2016) fail in the non-linear regime.

The HHG resulting from microscopic single-cell calculations using either TDDFT or SBEs (Fig. 2) displays a noisy spectrum with strong spectral contributions at all energies above the band-gap energy ($\epsilon_{\mathrm{gap}}^{\mathrm{LDA}}\approx 5.5$ eV $\approx 3.55\,\omega_{\mathrm{IR}}$). The lack of clear signatures of discrete harmonics in single-cell HHG spectra of an extended periodic system can be easily understood within a semiclassical picture Vampa et al. (2015b). Within the reduced real-space zone scheme of a single cell the electron and hole wave packets driven by the pulse traverse the cell many times meeting each other at a multitude of different recombination times and with different recombination energies. Only for high frequencies near the cut-off harmonic peaks become more clearly visible as the number of contributing recombination times is strongly reduced.

The broad quasi-continuous spectrum in the plateau region (black lines in Fig. 2) is at variance with a large number of solid-state HHG experiments finding a clean harmonic spectrum Ghimire et al. (2011); Schubert et al. (2014); Luu et al. (2015); Ndabashimiye et al. (2016). A recently employed phenomenological approach to “purify” the spectrum has invoked including very short dephasing times of the order of a fraction of an optical cycle ($T_{2}\approx T_{0}/4$ Vampa et al. (2014)) or of about 1 fs Luu et al. (2015); Luu and Wörner (2016) into the microscopic description. Such short $T_{2}$ raise, however, conceptual questions as to the origin of such ultrafast relaxation channels for electronic excitations.

Starting point of our quantitative estimate for $T_{2}$ are experimental data on the optical conductivity of dielectrics. We determine the frequency dependent conductivity $\sigma(\omega)$ from the experimental complex refractive index $\sqrt{\varepsilon(\omega)}=n(\omega)+i\kappa(\omega)$ Lide (2004) via

[TABLE]

For an impulsive broad-band excitation spectrum $F(t)\propto\delta(t)$ the time-dependent induced current density follows as

[TABLE]

using Ohm’s law. From Eq. (5) we deduce the frequency dependent decay constant $T_{2}(\omega)$ from fitting exponentials $\propto e^{-t/T_{2}(\omega_{i})}$ to the Gabor transform of Eq. (5) at various photon energies $\omega_{i}$ (Fig. 3).

For frequencies above the band gap (between 8 and 17 eV) we find $T_{2}(\omega_{i})\gtrsim 10$ fs considerably longer than previously assumed. It should be noted that this estimate is strictly valid only in the linear response regime. However, the present ab-initio multi-scale simulation allows to extend this estimate into the non-linear regime. The non-linear extinction coefficient at the IR driving frequency $\kappa(\omega_{\text{IR}},I_{0})$ can be directly determined from the attenuation of the driving field self-consistently propagated through the crystal. The resulting $\kappa$ is found to be strongly dependent on $T_{2}$ used in the simulation. Only for large $T_{2}$ ($>10$ fs) corresponding to small spectral broadening we find agreement with first experimental data Zimin et al. (2017) on non-linear $\kappa(\omega_{\text{IR}},I_{0})$ for diamond in the $10^{12}$ to $10^{13}$ W/cm2 regime. Using $T_{2}=10$ fs as lower bound in the SBEs (Fig. 2d and SI) shows that within a microscopic single-cell calculation, the lack of pronounced high harmonics in the plateau regime persists. Therefore, inclusion of mesoscopic effects of the extended system within a multiscale treatment is crucial to describe the build-up of a clean HHG spectrum.

The simulation including propagation (Eq. (3)) displays clearly visible harmonic peaks and a strongly reduced background even when dephasing is switched off ($T_{2}\rightarrow\infty$ in Fig. 4a,c, layer thickness $1\mu$m). Obviously, destructive interference between a multitude of electron-hole recombination events along the propagation path at different recombination times is key to the formation of a clean harmonic spectrum with peaks at odd harmonics. It is important to note that the peak positions agree with odd multiples of the blue-shifted transmitted driving frequency $\omega_{\text{t}}$ rather than the incident frequency $\omega_{\mathrm{IR}}$ (see Fig. 4c). This predicted non-linear blue shift $\delta\omega=\omega_{\text{t}}-\omega_{\mathrm{IR}}$ in diamond is another signature of the strong non-linear response of the solid and has been, indeed, experimentally observed in different materials Kuo and Ren (1992); Ghimire et al. (2011).

For the 1 $\mu$m thick crystal the directly transmitted pulse is followed by a second less intensive pulse from double internal reflection at the back and front boundaries of the slab. Due to its reduced intensity, the reflected pulse only weakly affects the total transmitted spectrum. The time delay between direct and reflected pulses ($\sim 17$ fs) allows for a direct and independent determination of the effective index of refraction of $n_{\mathrm{eff}}\approx 2.4$ (the real part of $\sqrt{\varepsilon(\omega)}$) in agreement with the experimental value at the carrier wavelength $\lambda=800$ nm Lide (2004).

Mesoscopic-scale effects include not only the propagation along the propagation direction but also the inhomogeneous field distribution within the focal spot in the transverse direction. The spot size is, typically, much smaller than the sample. Consequently, the intensity distribution of the laser pulse in the material has to be modeled adequately. The latter effect is included in this work by an average of the harmonic emission over a Gaussian profile $I(\rho)=I_{\text{max}}\exp(-\rho^{2}/2R^{2})$ of the driving laser pulse. To assess this effect independently, we calculate the resulting averaged current

[TABLE]

by a weighted sampling of the beam cross-section for an ensemble of 20 different intensities. Similar to the destructive interference due to propagation effects, also intensity averaging leads to spectra displaying more pronounced harmonic peaks with maxima at odd multiples of the driving frequency (Fig. S2 in SI).

Combining now the effects of destructive interference and, thus, dephasing due to propagation and transverse intensity averaging with dephasing due to microscopic decoherence processes with dephasing time $T_{2}=10$ fs (Fig. 4b and d) yield a well-defined harmonic spectrum in accord with experimental data. From the Gabor transform it is obvious that the effect of the finite microscopic dephasing time $T_{2}$ is primarily the damping of the (induced) post-pulse current $\vec{J}(t)$ on the time scale of a few tens of femtoseconds. It is, however, of minor importance for the formation of a clean harmonic spectrum.

In conclusion, we have presented the first ab-initio multi-scale simulation for solid-state high harmonics self-consistently combining the microscopic non-linear response in 3 dimensions within the framework of TDDFT or multiband Bloch equations with mesoscopic propagation and source distribution effects. Irrespective of the level of the underlying description, the microscopic simulation of the non-linear response of a single cell of the periodic system of diamond alone fails to yield a well-defined harmonic spectrum. It is the spatio-temporal distribution of the emission events on the mesoscopic scale that leads to the formation of a clean high-harmonic spectrum with pronounced peaks at odd harmonics of the blue-shifted driving frequency in the medium. This process can be viewed as the solid state analog to the shaping of harmonics by propagation in extended gaseous targets due to phase matching and suppression of the contributions from long trajectories Salières et al. (1998). Propagation and field inhomogeneity effects play, however, a much more prominent role at solid-state densities as in their absence the well-defined harmonic spectrum is replaced by a noisy quasi-continuum. Ultrafast microscopic dephasing rates of the order of $T_{2}\approx 1$ fs previously invoked are neither necessary nor justified for forming a well-defined harmonic spectrum. Instead, the present multi-scale simulation strongly suggests that the dephasing time should be at least one order of magnitude larger to yield the experimental optical conductivity at low intensities and the extinction coefficient in the non-linear regime of higher intensities.

The present work addressing diamond as a prototypical bulk dielectric has also wider important implications for other materials. The strong influence of propagation and field inhomogeneity effects in the dense medium on the non-linear optical response observed here is expected to be present irrespective of the details of lattice or electronic structure of the material. The blue shift of the harmonic spectrum found in the present multi-scale simulation for diamond has been, in fact, experimentally previously observed for other materials. The present findings may also contribute to disentangling of the relative importance of intraband and interband contributions to the harmonic generation. Since interband polarization is much more effectively suppressed by previously proposed ultrashort dephasing times than coherent intraband dynamics, the contribution of interband polarization may be significantly higher than expected when realistic dephasing times are employed. The present findings have also important implications for other ultra-fast processes, for example, probing of biomaterials by transmission of a broadband pulse Brinks et al. (2014). Extraction of information on the presence of specific molecules will require the disentangling of the non-linear molecular response from mesoscopic light transport effects.