Quantum beats in the polarization response of a dielectric to intense few-cycle laser pulses

Multiphoton strong-field excitation of valence-band electrons has been intensively investigated over the last several decades. Until recently, the main motivation for experimental research has been laser damage of optical materials [1–3]. Detailed time-resolved measurements of electron dynamics used to be out of reach, as the characteristic times of the most essential processes belong to the attosecond domain. The situation has recently changed due to the progress in the generation of few-cycle laser pulses and the development of attosecond science [4, 5]. In this context, several recent papers should be pointed out. Gertsvolf et al [6] observed that an elliptically polarized 40 fs pulse of near-infrared (NIR) radiation changes its ellipticity as it propagates through a few micrometres of fused silica. Based on their numerical calculations, these authors concluded that the change of ellipticity is due to sub-cycle, attosecond-scale excitation dynamics of valence-band electrons. Their conclusions were confirmed by Mitrofanov et al [7], who used a non-collinear pump–probe technique to observe an attosecond-scale modulation of the electron density created by a 5 fs NIR pulse in a SiO₂ sample. A different optical strong-field phenomenon was observed by Ghimire et al [8] in ZnO exposed to high-power few-cycle mid-infrared laser pulses: they observed the generation of high harmonics with properties that cannot be explained by the conventional nonlinear optics. The origin of these harmonics was linked to a combination of anharmonic electron motion and multiple Bragg reflections at the boundaries of the first Brillouin zone, also known as 'Bloch oscillations' [9]. Another manifestation of Bloch oscillations observed by the same team was a redshift in the optical absorption in ZnO [10]. Most recently, it was demonstrated that an intense few-cycle laser pulse not only excites charge carriers in SiO₂, but it also drives measurable currents that can be steered by controlling the carrier-envelope phase (CEP) [11] of the laser pulse [12]. The theoretical interpretation of that experiment was given in terms of the Wannier–Stark localization [13, 14] that can be viewed as a frequency-domain description of Bloch oscillations. These developments created a new perspective for studying electron dynamics in solids exposed to strong fields, where the emphasis is placed on attosecond-scale dynamics.

The foundations for the theoretical description of extremely nonlinear interaction with strong fields, where conventional perturbation theory breaks down, were laid by Keldysh in his seminal papers [15, 16] where multiphoton and tunnelling excitations [17], which are two extreme regimes of interband transitions, were distinguished by the Keldysh parameter

$$\begin{equation} \gamma = \frac{\omega_{\mathrm{L}} \sqrt{m E_{\mathrm{g}}}}{e F_0}. \end{equation} \tag{ 1 }$$

This parameter was introduced for a monochromatic laser field in the approximation of parabolic bands. In (1), ω_L is the laser frequency, m = (1/m_e + 1/m_h)⁻¹ is the reduced mass of an electron and a hole, E_g = ℏω_g is the bandgap, e is the absolute value of the electron charge and F₀ is the amplitude of the laser field in the medium. In the limit γ ≫ 1, where the field is weak and/or its frequency is large, an electron is excited from the valence band by absorbing an integer number of photons in a way that is well described by conventional perturbation theory. In the opposite extreme γ ≪ 1, where the field is strong and its frequency is sufficiently small, it is quantum tunnelling that promotes electrons to the conduction band. Obviously, this is the case if the external field is constant (ω_L = 0), and interband tunnelling as a concept was first introduced by Zener [18] for time-independent fields. In the tunnelling regime, the excitation rate is often assumed to be determined by the instantaneous value of the electric field F_L(t) [7].

Electron dynamics in the intermediate regime γ ∼ 1 are particularly complex—a laser field cannot be treated as a small perturbation in this regime and, at the same time, interband transitions cannot be described by pure tunnelling. A general expression for the cycle-averaged excitation rate was given by Keldysh in [16], but a completely satisfactory description of electron dynamics within an optical cycle of an ultrashort pulse is still missing. The dynamics can be modelled by solving the time-dependent Schrödinger equation (TDSE) in a suitable approximation, but there is no unambiguous way to distinguish between valence- and conduction-band electrons as long as the external field is present. This is analogous to the problem of determining the ionization rate of atoms and molecules in strong fields [19], for which many practical solutions exist, but there is no mathematically rigorous solution in the case of a few-cycle pulse. This motivated us to study the polarization response of a dielectric in the regime γ ∼ 1, as interband transitions directly affect polarization and, at the same time, polarization is a quantum-mechanical observable that determines all optical properties of a solid.

We are mainly interested in the strong-field regime where perturbation theory shows its shortcomings. In particular, the polarization response of a solid to a laser pulse is no longer described by a set of linear and nonlinear susceptibilities; therefore, quantum-mechanical simulations are necessary to model it. We solve the TDSE numerically in the velocity gauge, using the method described in [20]. For each value of the quasimomentum k, the electron wave function is represented in the basis of Bloch states $\left | \phi _{\boldsymbol {k}}^n\right \rangle$, which are evaluated by solving the single-electron stationary Schrödinger equation with an unperturbed Hamiltonian $\hat {H}_0$:

$$\begin{equation} \skew3\hat{H}_0\left| \phi_{\boldsymbol{k}}^n\right\rangle = E_{\boldsymbol{k}}^n \left|\phi_{\boldsymbol{k}}^n\right\rangle. \end{equation} \tag{ 2 }$$

Here, n is a band index. This equation is solved in the basis of N_max plane waves:

$$\begin{equation} \phi_{\boldsymbol{k}}^n({\boldsymbol{r}}) = \sum_{j=1}^{N_{\mathrm{max}}} C_{{\boldsymbol{k}},j}^n \exp[{\rm i} ({\boldsymbol{k}}+{\boldsymbol{K}}_j) {\boldsymbol{r}}], \end{equation} \tag{ 3 }$$

where K_j denotes reciprocal lattice vectors, each of which can be expressed as a combination of the primitive vectors {b₁,b₂,b₃}:

$$\begin{equation} {\boldsymbol{K}} = \sum_{\alpha=1}^3 m_\alpha {\boldsymbol{b}}_\alpha,\quad m_\alpha \in \mathbb{Z}. \end{equation} \tag{ 4 }$$

Ansatz (3) ensures that all Bloch states satisfy $\phi _{\boldsymbol {k}}^n({\boldsymbol {r}}+{\boldsymbol {R}}) = \phi _{\boldsymbol {k}}^n({\boldsymbol {r}}) \exp [{\mathrm { i}} {\boldsymbol {k}} {\boldsymbol {R}}]$, where R is a vector in the direct Bravais lattice. This property immediately follows from $\exp [{\mathrm { i}} {\boldsymbol {K}} {\boldsymbol {R}}] = 1$.

The coefficients in expansion (3), as well as the energies Eⁿ_k, are determined by diagonalizing the Hamiltonian $\hat {H}_0$ in the basis of plane waves $\exp [{\mathrm { i}} {\boldsymbol {K}}_j {\boldsymbol {r}}]$. We assume that the interaction of a particular electron with lattice ions and other electrons is described by a pseudopotential U(r) and expand this pseudopotential in the plane-wave basis:

$$\begin{equation} U({\boldsymbol{r}}) = \sum_{j=1}^{N_{\mathrm{max}}} U_j {\rm e}^{{\rm i} {\boldsymbol{K}}_j {\boldsymbol{r}}}, \end{equation} \tag{ 5 }$$

$$\begin{equation} U_j = \frac{1}{\Omega} \int_{\Omega} {\rm e}^{-{\rm i} {\boldsymbol{K}}_j {\boldsymbol{r}}} U({\boldsymbol{r}}) \,{\rm d}{\boldsymbol{r}}, \end{equation} \noindent \tag{ 6 }$$

where Ω is the volume of a unit cell. Writing the unperturbed Hamiltonian as

$$\begin{equation} \hat{H}_0 = -\frac{\hbar^2}{2 m_0} \nabla^2 + U({\boldsymbol{r}}) \end{equation} \noindent \tag{ 7 }$$

and substituting ansatz (3) into (2), we obtain the following eigensystem for the expansion coefficients Cⁿ_k,i and energies Eⁿ_k:

$$\begin{equation} \sum_{j=1}^{N_{\mathrm{max}}} \left( \frac{\hbar^2}{2 m_0} ({\boldsymbol{k}}+{\boldsymbol{K}}_j)^2 \delta_{i j} + U_{{\boldsymbol{K}}_i-{\boldsymbol{K}}_j} \right) C_{{\boldsymbol{k}},j}^n = E_{\boldsymbol{k}}^n C_{{\boldsymbol{k}},i}^n. \end{equation} \noindent \tag{ 8 }$$

Here, m₀ is the electron rest mass, and the eigensystem can be solved independently for each quasimomentum k. The solutions of this eigensystem are normalized to satisfy

$$\begin{equation} \left\langle \phi_{\boldsymbol{k}}^m | \phi_{\boldsymbol{k}}^n\right\rangle = \frac{1}{\Omega} \int_{\Omega} \left[\phi_{\boldsymbol{k}}^m({\boldsymbol{r}})\right]^* \phi_{\boldsymbol{k}}^n({\boldsymbol{r}})\,{\rm d}^3{\boldsymbol{r}} = \delta_{m n}, \end{equation} \noindent \tag{ 9 }$$

which implies that $\sum _j |C_{{\boldsymbol {k}},j}^n|^2 = 1$. Here and in the following, a scalar product $\left \langle \cdot | \cdot \right \rangle$ implies integration over one unit cell.

Having evaluated Bloch states $\left | \phi _{\boldsymbol {k}}^n\right \rangle$, we use them as a basis to solve the TDSE and thus model the interaction of electrons with a laser pulse:

$$\begin{equation} {\rm i} \hbar \frac{\partial\,}{\partial t} \left| \psi_{\boldsymbol{k}}(t)\right\rangle = (\hat{H}_0 + \hat{H}_{\mathrm{int}}(t)) | \psi_{\boldsymbol{k}}(t)\rangle, \end{equation} \tag{ 10 }$$

where $\hat {H}_{\mathrm {int}}(t)$ is the interaction Hamiltonian. We perform our simulations in the velocity gauge, where an external electric field is described by the vector potential A(r,t). We also employ the dipole approximation: A(r,t) ≡ A(t). In SI units, this corresponds to using the following form of the interaction Hamiltonian:

$$\begin{equation} \hat{H}_{\mathrm{int}}(t) = \frac{e}{m_0} {\boldsymbol{A}}(t) \hat{{\boldsymbol{p}}}. \end{equation} \noindent \tag{ 11 }$$

Note that in CGS units the right-hand side of (11) must be divided by the vacuum speed of light c. The substitution of the ansatz

$$\begin{equation} \left| \psi_{\boldsymbol{k}}(t)\right\rangle = \sum_{n=1}^{N_{\mathrm{max}}} \alpha_{\boldsymbol{k}}^n(t) \left| \phi_{\boldsymbol{k}}^n\right\rangle \end{equation} \noindent \tag{ 12 }$$

into (10) leads to the following system of coupled differential equations:

$$\begin{equation} {\rm i} \hbar \frac{\partial\,\alpha_{\boldsymbol{k}}^q(t)}{\partial t} = E_{\boldsymbol{k}}^q \alpha_{\boldsymbol{k}}^q(t) + \frac{e}{m_0} {\boldsymbol{A}}(t) \sum_{l=1}^{N_{\mathrm{max}}} {\boldsymbol{p}}_{\boldsymbol{k}}^{q l} \alpha_{\boldsymbol{k}}^l(t). \end{equation} \tag{ 13 }$$

Here,

$$\begin{equation} {\boldsymbol{p}}_{\boldsymbol{k}}^{q l} = \left\langle \phi_{\boldsymbol{k}}^q\left| \hat{{\boldsymbol{p}}} \right| \phi_{\boldsymbol{k}}^l\right\rangle = \hbar \sum_{j=1}^{N_{\mathrm{max}}} ({\boldsymbol{k}}+{\boldsymbol{K}}_j) \left(C_{{\boldsymbol{k}},j}^q\right)^* C_{{\boldsymbol{k}},j}^l \end{equation} \tag{ 14 }$$

are the matrix elements of the momentum operator $\hat {{\boldsymbol {p}}}$. An important advantage of the velocity gauge is that, as long as the dipole approximation holds and electron scattering is neglected, there is no coupling between different values of k. Mathematically, this is a consequence of

$$\begin{eqnarray*} &&\int_{\mathbb{R}^3} \left(\phi_{{\boldsymbol{k}}'}^q({\boldsymbol{r}})\right)^* \hat{{\boldsymbol{p}}}\, \phi_{\boldsymbol{k}}^l({\boldsymbol{r}}) \,{\rm d}^3{\boldsymbol{r}} = \left\langle \phi_{{\boldsymbol{k}}'}^q\left| \hat{{\boldsymbol{p}}} \right| \phi_{{\boldsymbol{k}}}^l\right\rangle \sum_{j} {\rm e}^{{\rm i} ({\boldsymbol{k}}-{\boldsymbol{k}}') {\boldsymbol{R}}_j } \propto \delta_{{\boldsymbol{k}} {\boldsymbol{k}}'}, \end{eqnarray*}$$

where the summation is performed over all the vectors of the Bravais lattice R_j, and we assume that both k and k' belong to the first Brillouin zone. Due to this property, equations (13) can be solved independently for each quasimomentum k. According to our experience, the downside of the velocity gauge is a larger number of bands required for convergence, as compared to length-gauge simulations, where different values of k are coupled to each other [21, 22], but even the two-band approximation may be sufficient [23].

At the beginning of a simulation, all electrons are supposed to reside in the valence bands (VB), while all the conduction bands are empty. We solve equations (13) for each of the electrons independently and then add single-electron contributions together to evaluate physical observables such as current density or polarization. Obviously, our approach is only applicable if (i) lattice vibrations have no significant effect on the investigated dynamics, (ii) effects related to electron correlation and electron–hole interaction are not important and (iii) distortions of the pseudopotential U(r) upon electronic excitations are negligible. While these approximations obviously have their limitations, single-electron models proved to be useful in studying basic physics related to phenomena occurring on the attosecond and few-femtosecond time scales [20, 24, 25].

The initial condition for solving equations (13) for a particular k is

$$\begin{equation} \alpha_{\boldsymbol{k}}^{q}(t_{\mathrm{min}}) = \delta_{q n}. \end{equation} \noindent \tag{ 15 }$$

Here, n is a band where the electron was before the interaction with a laser pulse.

In this paper, we are most interested in the macroscopic electric current density j(t) induced by the laser pulse

$$\begin{equation} {\boldsymbol{j}}(t) = \sum_{n \in \mathrm{VB}} \int {\boldsymbol{j}}_{{\boldsymbol{k}},n}(t)\,{\rm d}^3{\boldsymbol{k}}, \end{equation} \tag{ 16 }$$

where the integral is taken over the first Brillouin zone and j_k,n(t) represents single-electron contributions to the current density averaged over a unit cell:

$$\begin{equation} \fl {\boldsymbol{j}}_{{\boldsymbol{k}},n}(t) = -\frac{e}{m_0} \frac{1}{\Omega} \int_{\Omega} {\rm d}^3{\boldsymbol{r}}\, \left\{ {\rm Re}\left[\psi_{{\boldsymbol{k}},n}^*({\boldsymbol{r}},t) \hat{{\boldsymbol{p}}} \psi_{{\boldsymbol{k}},n}({\boldsymbol{r}},t)\right] +e {\boldsymbol{A}}(t) |\psi_{{\boldsymbol{k}},n}({\boldsymbol{r}},t)|^2 \right\}. \end{equation} \noindent \tag{ 17 }$$

Here and in the following, we add index n to quantities related to single-electron wave functions in order to indicate the initial band. Thus, j_k,n(t) represents the contribution from all the bands to j(t) in a calculation where the electron initially occupied band n. The term containing the vector potential is required in the velocity gauge.

With the aid of (12), the single-electron current density can be expressed via probability amplitudes α^q_k(t):

$$\begin{equation} {\boldsymbol{j}}_{{\boldsymbol{k}},n}(t) = -\frac{e}{m_0} \left({\rm Re}\left[ \sum_{q,l} \left( \alpha_{{\boldsymbol{k}},n}^q(t) \right)^* \alpha_{{\boldsymbol{k}},n}^l(t) {\boldsymbol{p}}_{\boldsymbol{k}}^{q l} \right] + e {\boldsymbol{A}}(t) \right). \end{equation} \tag{ 18 }$$

Three different physical effects determine the induced currents: polarization due to the displacement of valence-band electrons, the light-driven motion of conduction-band electrons and quantum beats due to coherent superpositions of valence- and conduction-band states. In the presence of a strong external electric field, these contributions cannot be fully separated from each other, but the overall induced current density j(t) is defined unambiguously. It is this induced current that determines the change of the refractive index and the generation of new frequency components. However, it is conventional to describe optical response in terms of polarization, which we define as

$$\begin{equation} {\boldsymbol{P}}(t) = \int_{-\infty}^t {\boldsymbol{j}}(t')\,{\rm d} t'. \end{equation} \tag{ 19 }$$

Note that the polarization defined this way is due to both bound and free electrons. Ordinarily, P(t) would be assigned to bound electrons, while j(t) would describe the motion of free electrons. We make no attempt to make this distinction because, to the best of our knowledge, there is no rigorous way to distinguish between bound and free electrons as long as the external field is present [19]. It is only after the laser pulse that projecting a wave function onto conduction-band Bloch states yields physically relevant excitation probabilities.

For simplicity, we perform our numerical simulations in one spatial dimension x, which we assume to be the polarization direction of the laser field. Most of the equations in this section can be transformed to their one-dimensional (1D) forms by replacing all vector quantities and 3D integrals with scalar quantities and 1D integrals: r → x, k → k. The only exception is (16), where this procedure would result in j(t) being measured in wrong units. Reducing the dimensionality of our original problem, we essentially assume that physical observables do not depend on the coordinates y and z, so that the integral over the first Brillouin zone in (16) reduces to

$$\begin{equation} j(t) = \eta \sum_{n \in \mathrm{VB}} \int j_{k,n}(t)\,{\rm d} k, \end{equation} \tag{ 20 }$$

where the factor η, measured in m⁻², accounts for the integration over k_y and k_z. In practice, the value of η should be chosen to approximate some known optical properties of the solid, such as its refractive index or absorption.

From this point on, we will write our equations in atomic units (au), unless stated otherwise. In atomic units, ℏ = e = m₀ = 1, a unit of energy is equal to 27.21 eV, and a unit of electric field is equal to 5.142 × 10¹¹ V m⁻¹. Because ℏ = 1, we will interchangeably use the words 'energy' and 'frequency'.

We adjusted our model potential U(x) to obtain a bandgap close to that of SiO₂: 8.95 eV [26]. To this end, we set the lattice constant to a₀ = 0.5 nm = 9.45 au, which is one of the lattice constants in α-quartz, and use the following expression for the central unit cell:

$$\begin{equation} U(x) = -0.7 \left(1 + \tanh[x+0.8]\right) \left(1 + \tanh[-x+0.8]\right),\quad |x| \leqslant a_0/2. \end{equation} \noindent \tag{ 21 }$$

Outside the central unit cell, the potential is continued periodically: U(x + a₀) = U(x). The band structure that corresponds to this potential is shown in figure 1(a). At k = 0 the energies Eⁿ_k of the first four bands are −35.48, −8.33, 0.62 and 13.13 eV. We regard the first two bands with negative values of Eⁿ_k as valence bands (the lowest valence band is not shown in figure 1(a)), while all the bands with Eⁿ_k > 0 are regarded as conduction bands.

Zoom In Zoom Out Reset image size

We used the following expression for the vector potential of the (linearly polarized) laser field acting on electrons in the solid:

$$\begin{equation} A(t) = -\theta(\tau_{\mathrm{L}}-|t|) \frac{F_0}{\omega_{\mathrm{L}}} \cos^4\left(\frac{\pi t}{2 \tau_{\mathrm{L}}}\right) \sin(\omega_{\mathrm{L}} t + \varphi_{\mathrm{CEP}}). \end{equation} \tag{ 22 }$$

Here, F₀ is the amplitude of the electric field, ω_L is the central frequency of the laser pulse, φ_CEP is its CEP, θ(x) is the Heaviside step function and τ_L is related to the full-width at half-maximum (FWHM) of |A(t)|² by ${\mathrm { FWHM}} = 4 \arccos (2^{-1/8}) \tau _{\mathrm {L}} / \pi \approx 0.5224 \tau _{\mathrm {L}}$. The FWHM of |A(t)|² is very close, although not precisely equal, to the FWHM of |F(t)|². The external electric field acting on electrons is, by definition,

$$\begin{equation} F(t) = -A'(t). \end{equation} \tag{ 23 }$$

Note that, in this work, we do not make any attempt to evaluate the screening field created by the displacement of electrons—this field is assumed to be a part of A(t) and F(t).

For a weak pulse, the polarization response is linear, so that, in the frequency domain, $\tilde {P}(\omega ) = \chi (\omega ) \tilde {F}(\omega )$, and the refractive index can be evaluated from the linear susceptibility χ(ω):

$$\begin{equation} n(\omega) = \sqrt{1 + \chi(\omega)}. \end{equation} \tag{ 24 }$$

In figure 1(b), we compare the refractive index of fused silica with n(ω) evaluated using our model. The refractive index is plotted against the laser wavelength λ = 2πc/ω. For these simulations, we used several pulses with F₀ = 10⁻⁵ au = 5 × 10⁶ V m⁻¹, FWHM = 4π/ω_L, and values of ω_L that allowed us to cover the spectral range presented in figure 1(b). Normalizing the polarization response, we set η = 0.111 au in (20) and use this value henceforth. Given the simplicity of our model, we find that the agreement with the measured refractive index is very satisfactory.

The next section reports on simulations with much more intense fields, where multiphoton excitations populate conduction bands. For those simulations, we had to use 15 bands in order to achieve numerical convergence in solving the TDSE. A smaller number of bands results in discrepancies that are most visible in the spectral range occupied by low-order harmonics, while increasing the number of bands to 20 has a negligible effect on the polarization response even for the most intense pulses that we used in our modelling.

For all the simulations in this section, we used pulses with the FWHM equal to three optical oscillations: FWHM = 3 × 2π/ω_L. The electric field F(t) of such a pulse is shown in figure 2(a) by the dotted blue line plotted against the number of optical cycles ω_Lt/(2π). The solid red line in this figure shows the polarization response P(t) evaluated according to (19) for the case where the central frequency is equal to a quarter of the bandgap (ω_L = ω_g/4, λ_L = 2πc/ω_L = 557 nm), and the peak intensity is equal to I_L = 3 × 10¹³ W cm⁻², which corresponds to an amplitude of the electric field of F₀ = 0.029 au = 1.5 V Å⁻¹. The scales chosen for the figure emphasize that, at the leading edge of the laser pulse, P(t) is proportional to F(t). At later times, the polarization response becomes increasingly nonlinear and high-frequency oscillations appear, which persist even after the laser pulse ends. The frequency of these oscillations is close to the bandgap frequency ω_g = E_g/ℏ. Thus, these are quantum beats appearing due to the presence of coherent superpositions of valence- and conduction-band states. This part of the polarization response will be the main topic of our further discussion. In order to distinguish it from conventional harmonics due to χ⁽³⁾, χ⁽⁵⁾, χ⁽⁷⁾,... nonlinearities, we chose the central laser frequency to be an even divisor of the bandgap, taking advantage of the fact that harmonics of even orders are absent if the potential U(x) possesses the inversion symmetry U(−x) = U(x), like our model potential (21). Note that, due to the same symmetry, the four-photon transition from the uppermost valence band to the lowest conduction band at k = 0 is forbidden, while interband tunnelling can efficiently populate conduction-band states in the middle of the Brillouin zone. Indirectly, this is confirmed by the fact that the fast oscillations in figure 2(a) only appear at very high intensities. At a peak intensity of I_L = 10¹³ W cm⁻², which is just three times smaller than the one used for figure 2, the amplitude of quantum-beat oscillations is smaller by approximately a factor of 400, which corresponds to a drop of their intensity by five orders of magnitude.

Zoom In Zoom Out Reset image size

We investigate the quantum-beat part of the polarization response by means of the wavelet analysis using the Morlet wavelet [28]:

$$\begin{eqnarray} &&S(t,\omega) = \omega \left| \int_{-\infty}^\infty P(t') W\bigl(\omega (t' - t)\bigr)\,{\rm d} t' \right|^2, \end{eqnarray} \noindent \tag{ 25 }$$

$$\begin{eqnarray} &&W(x) = \frac{1}{\sqrt{\tau_W}} \exp\left[{\rm i} x - \frac{x^2}{2 \tau_W^2}\right]. \end{eqnarray} \noindent \tag{ 26 }$$

The parameter τ_W determines the size of the temporal window, and it must be chosen to provide an optimal compromise between the temporal and spectral resolution. Our choice was τ_W  = 3ω_g/ω_L.

The outcome of this analysis, S(t,ω), is shown as a false-colour diagram in figure 2(b), the bandgap energy being indicated by a dashed white line. For these parameters, the quantum-beat signal completely dominates the third-order harmonic, which is barely visible in the plot. In this time–frequency analysis, the third-harmonic signal is centred at 6.7 eV, and it is temporally confined to the FWHM of the light pulse: ω_L|t|/(2π) ≲ 1.5.

One of the most striking features observed in figure 2(b) is the fact that the quantum-beat signal initially appears at frequencies exceeding ω_g by roughly 2 eV/ℏ. The instantaneous frequency of quantum beats then gradually decreases, approaching ω_g at the tailing edge of the laser pulse. In other words, the quantum-beat signal is negatively chirped—its lower-frequency components are delayed with respect to the higher-frequency ones. This cannot be explained by the presence of the fifth-order harmonic—even though the fifth harmonic occupies the relevant spectral range (5ℏω_L = 11.2 eV), it must be confined to an even smaller temporal range than the third harmonic, so that S(t,ω) can be considered free from its contribution for ω_Lt/(2π) ≳ 1, where most of the quantum-beat signal is observed.

We seek an explanation for the chirp of the quantum-beat signal in the laser-driven motion of charge carriers. The average kinetic energy of an electron–hole pair is the ponderomotive energy. For a monochromatic field with an amplitude F₀, the ponderomotive energy is given, in the approximation of parabolic bands, by

$$\begin{equation} U_{\mathrm{p}} = \frac{F_0^2}{4\,m \omega_{\mathrm{L}}^2}. \end{equation} \tag{ 27 }$$

Here, m is the reduced mass of an electron and a hole. We obtain the effective masses of electrons and holes by fitting parabolas to the lowest conduction and uppermost valence bands in the region |k| < π/(2a₀), which yields m_e = 0.34 au and m_h = 3.61 au. These values are in good agreement with those that can be found in the literature [29, 30]. Thus, the reduced effective mass m in (27) is equal to m = (1/m_e + 1/m_h)⁻¹ = 0.35 au. The cycle-averaged total energy of an electron–hole pair is equal to E_g + U_p(t), where U_p(t) is to be evaluated by replacing F₀ in (27) with the pulse envelope. This quantity is plotted as a red line in figure 2(c). It is fairly close to the energies where S(t,ω) has local maxima with respect to ω. We observed this kind of negative chirp in numerous simulations with different parameters of the laser pulse, some of which are shown in figure 3. Therefore, we conclude that the frequency of quantum beats exceeds ω_g by approximately the ponderomotive energy.

Zoom In Zoom Out Reset image size

Figure 2(d) presents a different kind of analysis that highlights the quantum-beat part of the polarization response. In order to suppress the third- and fifth-order harmonics, we multiplied the Fourier transform of P(t) with a soft spectral window that cuts all frequency components that are further than ω_L/2 from the bandgap frequency:

$$\begin{equation} \tilde{w}(\omega) = \theta\left(\frac{\omega_{\mathrm{L}}}{2} - \biggl||\omega|-\omega_{\mathrm{g}}\biggr|\right) \cos^2\left(\pi \biggl||\omega|-\omega_{\mathrm{g}}\biggr| / \omega_{\mathrm{L}}\right). \end{equation} \noindent \tag{ 28 }$$

The grey line in figure 2(d) shows the inverse Fourier transform of $\tilde {w}(\omega ) \tilde {P}(\omega )$. In other words, it represents the part of the polarization response that oscillates at frequencies close to the bandgap frequency. The envelope of these fast oscillations is shown as a solid red line, and the dotted blue line depicts the square of the laser field. From this analysis, it is obvious that the quantum-beat part of the polarization response experiences its most rapid increase at a time that lies about two optical cycles after the peak of the laser pulse. This is counter-intuitive, as the quantum-beat signal is due to interband excitations, and the rate of the interband excitations is expected to reach its peak at the peak of the laser pulse. Another prominent feature apparent from figure 2(d) is the decay of quantum beats that begins at the trailing edge of the laser pulse. As our model does not account for any relaxation phenomena (such as electron–phonon interaction or spontaneous emission of radiation), the origin of this decay lies in the dephasing of dipole oscillators associated with each pair of coherently populated valence- and conduction-band states. Each such coherent superposition creates a current j_k(t) described by (18). If the laser pulse ends at a time t₁, then $\alpha _k^q(t) = \alpha _k^q(t_1) \exp [-{\mathrm { i}} (t-t_1) E_k^q]$, so that

$$\begin{equation} j_k(t) = -{\rm Re}\left[ \sum_{q,l} \left( \alpha_k^q(t_1) \right)^* \alpha_k^l(t_1) p_k^{q l} {\rm e}^{{\rm i} (t-t_1) (E_k^q-E_k^l)} \right]. \end{equation} \noindent \tag{ 29 }$$

For brevity, we have omitted the index specifying the initial band of the electron. For a given pair of bands q ≠ l, (29) describes a current oscillating with a k-dependent frequency E^q_k − E^l_k. Because of this dependence, the net current density (20) attenuates as most currents j_k(t) get out of phase with each other.

Figure 3 presents the wavelet analysis of polarization response functions for several other values of the central laser frequency ω_L. The dashed white lines depict E_g + U_p(t). In figures 3(a) and (b), these lines closely follow local maxima of S(t,ω). For smaller laser frequencies, the ponderomotive potential U_p significantly exceeds the width of the lowest conduction band (see figure 1(a)). This is the regime where Bloch oscillations [31] must play an important role [8].

In figure 3(a), where valence-band electrons can be excited by five-photon absorption, we made an exception from our choice to use those values of ω_L that are even divisors of ω_g. Still, the time–frequency analysis looks very similar to that for ω_L = ω_g/4, shown in figure 2(b), and the contribution from the fifth harmonic generation is so weak that it is practically invisible in our colour scheme. For all the cases shown in figure 3, we observe the same qualitative features as those in figure 2(b): the frequency components of the polarization response that lie above the bandgap appear delayed with respect to the peak of the laser pulse, and the instantaneous frequency of these oscillations decreases with time, eventually approaching ω_g. We also observe a rapid decrease in the magnitude of S(t,ω).

To further understand the origin of these trends, we plot the sum of all conduction-band populations at the end of our simulations in figure 4. These are the probabilities of finding an electron with a given quasimomentum k in one of the conduction bands at a time t > τ_L where the laser field is absent (since we do not account for relaxation phenomena, populations do not change in the absence of external fields). The distribution of excited electrons dramatically changes as ω_L decreases from ω_g/4 to ω_g/10: the pronounced peak at k = 0 disappears, and many irregular peaks appear across the Brillouin zone. At the same time, the distributions become increasingly sensitive to the CEP of the laser pulse. We consider two physical phenomena that may be responsible for this trend: interband tunnelling and Bloch oscillations.

Zoom In Zoom Out Reset image size

The first phenomenon, interband tunnelling, represents an extreme regime of interband excitations where the Keldysh parameter γ is much smaller than 1. In this regime, interband transitions should be viewed not as a result of absorbing a certain number of photons, but as a result of quantum tunnelling [18]. Also, in the tunnelling regime, the CEP is an important parameter because it is the field of the laser pulse, rather than its envelope, that controls the tunnelling rate. According to figure 4 and all the other similar simulations we have performed, the onset of the CEP dependence coincides with γ decreasing below a value approximately equal to 0.7. Furthermore, the tunnelling picture provides an intuitive explanation for the multitude of peaks in figures 4(c) and (d): they may be attributed to the interference among electron wave packets launched to conduction bands within different half-cycles of the laser pulse, which is analogous to the above-threshold ionization spectra of atoms exposed to intense laser pulses [32].

What the tunnelling picture alone cannot explain is the disappearance of the peak at k = 0. Indeed, all analytical models for the tunnelling rate predict that it should rapidly decrease as the bandgap increases [16, 18, 33], and our model potential has the smallest bandgap at k = 0. Therefore, the motion of excited electrons in conduction bands must play a significant role [23], Bloch oscillations being the most prominent manifestation of this motion. An electron excited to the bottom of a band at a zero-crossing of the vector potential (A(t₀) = 0) may reach the boundary of the first Brillouin zone if, at some later moment t₁, the condition |A(t₁)| = π/a₀ is satisfied. Thus, Bloch oscillations should be considered important if the amplitude of the vector potential F₀/ω_L is comparable to π/a₀. Introducing the Bloch frequency

$$\begin{equation} \omega_{\mathrm{B}} = \frac{e a_0 F_0}{\hbar}\, {\rm (SI~units)}, \end{equation} \tag{ 30 }$$

this condition can be expressed as ω_B/ω_L = π. The labels in figures 4(a)–(d) provide values of ω_B/ω_L in units of π for all the four simulations. The appearance of the CEP dependence and the complex structures in the distribution of conduction-band electrons coincide not only with γ becoming smaller than 1, but also with ω_B/ω_L becoming larger than π. In fact, this is not a coincidence, as

$$\begin{equation} \frac{\omega_{\mathrm{B}}}{\omega_{\mathrm{L}}} = \frac{a_0 \sqrt{m E_{\mathrm{g}}}}{\hbar \gamma}\quad {\rm (SI~units)}, \end{equation} \noindent \tag{ 31 }$$

as follows from equations (1) and (30). For our parameters, ω_B/ω_L ≈ 1.02π/γ, so that the onset of tunnelling coincides with the onset of Bloch oscillations, irrespective of the laser pulse parameters.

Without interband excitations, Bloch oscillations obviously cannot occur, but is it the nature of interband transitions or the Bloch oscillations themselves that are mainly responsible for the observed CEP dependences? While we cannot answer this question with certainty, we can argue that interband tunnelling alone is not sufficient for a polarization response to be CEP dependent. Indeed, tunnelling must be the dominant excitation mechanism in the parameter regime of figure 4(a) because, in spite of the forbidden four-photon transitions at k = 0, the bottom of the conduction band is strongly populated after the laser pulse. However, the polarization response is almost CEP independent. Therefore, our simulations with ω_L = ω_g/4 provide an example where the dominance of interband tunnelling over (perturbative) multiphoton excitations does not lead to a significant CEP dependence.

The quantum-mechanical observables that we have discussed so far cannot be directly measured in a realistic experiment. The most accessible quantity for experiments is the spectrum of emitted radiation. In order to evaluate such spectra, one needs to account for absorption and phase-matching [9], which is beyond the scope of this work. Nevertheless, important conclusions about these spectra can be made by analysing the Fourier transform of the polarization response P(t) that we evaluate with our model. In figure 5, we plot polarization spectra for the same simulations as those presented in figure 4. We see that the spectra are also sensitive to the CEP of the laser pulse. Similarly to excitation probabilities, the CEP dependence is very weak for ω_L = ω_g/4, but it becomes very pronounced already at ω_L = ω_g/6, and it affects not only the quantum-beat part of the polarization response, but also the conventional harmonics of third, fifth, etc orders.

Zoom In Zoom Out Reset image size

Consistent with the results shown in figure 3, the quantum-beat signal rapidly decreases with decreasing ω_L. At the same time, according to figure 4, the excitation probabilities are comparable in all the four cases. Consequently, the decrease in the polarization response is due to dephasing. Most importantly, we used longer pulses for longer wavelengths in order to have the same number of optical cycles in all the pulses. As a result, dephasing played a more significant role for the 14 fs pulse with ω_L = ω_g/10 than for the 5.5 fs pulse with ω_L = ω_g/4. In addition to that, the dephasing time is inversely proportional to the spectral width of the excitation, and excited electrons occupy a larger part of the Brillouin zone in the tunnelling regime (see figure 4).

We have analysed the polarization response of a model dielectric resembling SiO₂ to few-cycle laser pulses that are strong enough to efficiently excite electrons from valence to conduction bands. These interband transitions create coherent superpositions of states that manifest themselves as quantum beats in the polarization. Obviously, the dynamics of interband excitations and dephasing determine the properties of the quantum-beat part of the polarization response. It is less obvious that the instantaneous frequency of quantum-beat oscillations changes with time, approaching the bandgap frequency as the laser field attenuates. We find that, during the laser pulse, the quantum-beat frequency (averaged over a laser cycle) is larger than the bandgap frequency by U_p/ℏ, unless the ponderomotive energy U_p exceeds the width of the lowest conduction band. Thus, the emitted burst of high-frequency radiation is negatively chirped. We have also observed that quantum beats appear in the polarization response significantly delayed with respect to the peak of a laser pulse, even though most electron–hole pairs must appear when the laser field is strongest. Apparently, the effect of a strong field on charge carriers suppresses single-photon transitions from conduction- to valence-band states. A satisfactory explanation of this phenomenon probably requires the development of an analytical theory, which may be triggered by this work.

For a dielectric such as SiO₂, the onset of interband tunnelling (γ ∼ 1) coincides with the onset of Bloch oscillations (ω_B ∼ πω_L). In the tunnelling regime, we have observed that the nonlinear polarization and the final distribution of charge carriers are sensitive to the CEP of a laser pulse. In particular, the spectrum of the emitted high-frequency radiation depends on the CEP, which should be possible to observe in experiments. Such a measurement would have to overcome obstacles related to the strong absorption above the bandgap and the phase mismatch below the bandgap [9]. Also, it must be mentioned that the coherent few-femtosecond pulse radiated by quantum beats spectrally overlaps with the incoherent radiation due to fluorescence, which our model does not account for. Some ultrashort temporal gating may therefore be necessary to isolate the radiation due to quantum beats. Still, in spite of all these obstacles, such measurements should be feasible at the current stage of technology.

The authors gratefully acknowledge fruitful discussions with Professor F Krausz and Professor K Yabana. This work was supported by the DFG Cluster of Excellence: Munich-Centre for Advanced Photonics.