Pulsed high harmonic generation of light due to pumped Bloch oscillations in noninteracting metals

# Pulsed high harmonic generation of light due to pumped Bloch oscillations in noninteracting metals

## Abstract

We derive a simple theory for high-order harmonic generation due to pumping a noninteracting metal with a large amplitude oscillating electric field. The model assumes that the radiated light field arises from the acceleration of electrons due to the time-varying current generated by the pump, and also assumes that the system has a constant density of photoexcited carriers, hence it ignores the dipole excitation between bands (which would create carriers in semiconductors). We examine the circumstances under which odd harmonic frequencies would be expected to dominate the spectrum of radiated light, and we also apply the model to real materials like ZnO, for which high-order harmonic generation has already been demonstrated in experiments.

Nonequilibrium, high harmonic generation, Bloch oscillations, zinc oxide
###### pacs:
78.47.je,78.20.-e,72.20.Ht

## I Introduction

The theory of Bloch oscillations has a long history, dating back to the original work by Bloch bloch () and Zener zener (). It remains one of the most difficult phenomena to observe in real materials due to the fast timescale for the oscillations, and it has been seen predominantly in semiconductor superlattices superlattice () and in ultracold atomic systems when placed on optical lattices optical_lattice (). We argue here that similar phenomena can be observed by examining the radiated light from solids that are pumped by high intensity laser pulses with femtosecond pulse widths kemper (); golde ().

We describe such laser pulses as spatially-uniform time-varying electric fields, using a semiclassical approach. In this sense, we ignore all magnetic field effects, and essentially violate Maxwell’s equations for the electromagnetic fields. This approach will be accurate if the magnetic field effects on the solid are small, which should be true for rapidly oscillating fields, because the electron spin will not be able to respond to the magnetic field on such short time scales. The large electric field will cause the electrons in the material to accelerate, giving rise to a time-dependent current. Assuming the current on the surface is proportional to the bulk current, we will see electromagnetic radiation coming off the material with a signal proportional to the time derivative of the current. We will use this to directly determine the spectra of the radiated light.

Experimental investigations of high harmonic generation in bulk crystals are being pursued, including recent work on irradiating ZnO crystals, where many tens of harmonics of the mid-infrared driving pulse were measured zno_exp (). For most crystallographic orientations, only the odd harmonics of the fundamental frequency appeared. We will see below that the occurence of only odd harmonics is to be expected for systems that satisfy time-reversal symmetry.

In the next section, we discuss the general formalism, and apply it to simple model systems, which can be thought of as tight-binding models for the bandstructure. In Section III, we use numerics to examine the high harmonic generation in simple models, and in Section IV, we use numerics to examine the high harmonic generation in ZnO and compare the numerical results with those of experiment. We conclude in the final section.

## Ii Formalism

We begin with a system described by a bandstructure, which we denote by . The bandstructure is a periodic function of the Brillouin zone. If we use a tight-binding description for the bands, then the bandstructure can generically be written as

 ϵ(k)=−∑δtδeik⋅δ, (1)

for -orbitals, where is a translation vector to a neighboring site (the sum need not stop at just nearest neighbors), and is the corresponding hopping integral. The band index, if there are multiple bands, has been suppressed. For models with higher angular momentum orbitals, the tight-binding form for the bandstructure is more complicated and requires diagonalization of a matrix whose size is equal to the number of different types of orbitals per unit cell (times two in the presence of magnetic or spin-orbit interactions).

We work in the Hamiltonian gauge, where the electric field is described by a spatially uniform vector potential , via

 E(t)=−1c∂A(t)∂t, (2)

with the speed of light. The electric field is then introduced into the bandstructure with the so-called Peierls substitution which takes , with the electron charge. This substitution properly describes the way in which electrons are accelerated in each band, but does not take into account the dipole coupling between different types of bands (like between and bands). Such a direct Zener tunneling term can be included, but will complicate the analysis, and we ignore it in the results we describe here, for simplicity. We expect that the results one finds if we include such terms will be similar, because the main effect of the laser pulse is to accelerate the electrons within a given band in metallic systems, and excitations between bands become more important for semiconductors or insulators. But even in that case, if we assume most of the photoexcitation occurs when the field amplitude is the highest, then the results of this work will be identical except for a factor of 2, which does not affect any of the results we present below.

The Hamiltonian of the system then becomes

 H(t)=∑kσϵ(k−eA(t)ℏc)c†kσckσ (3)

where () is the electron creation (annihilation) operator for an electron with momentum and spin (here we are assuming, for simplicity, a single band model and no spin-orbit coupling, so is a good quantum number). This is a time-dependent Hamiltonian, but it commutes with itself at different times, which makes its time evolution particularly easy to solve.

The system starts off in equilibrium at early times and the electric field pulse is centered around . Hence, the initial distribution of the electrons is given by the Fermi-Dirac distribution

 f(ϵ−μ)=11+exp[β(ϵ−μ)] (4)

in the distant past. Here, is the chemical potential, and is the inverse temperature. We work with an isolated system, so the time dynamics can be solved directly in the Heisenberg picture. The current satisfies

 J(t)=eℏ∑kσ∇kϵ(k−eA(t)ℏc)f(ϵk−μ), (5)

which simply sums the field-dependent velocity (), weighted by the initial distribution of the electrons and multiplied by the electric charge. This formula also follows from a Kadanoff-Baym-Keldysh formalism freericks_nonint () where one can exactly solve for the Green’s functions, evaluate them at equal time and sum the electron occupancy at weighted by the field-dependent velocity.

If we assume that the radiated light arises from the acceleration of the electrons, then the light spectrum is proportional to

 ∣∣∣∫dteiωtddtJ(t)∣∣∣2=∣∣∣ω∫dteiωtJ(t)∣∣∣2. (6)

If we have a pulse with a shape given by

 A(t)=−A0cos(ω0t+ϕ)exp(−t2/t20)t0√π, (7)

then we are guaranteed that the field has no dc component (required of optical pulses) and if is somewhat larger than the period of the oscillations , then the pumped electric field has a reasonably well defined fundamental frequency ( is an adjustable phase shift of the oscillations relative to the gaussian envelope). While we choose a gaussian form for the envelope function, any normalized slowly varying function which vanishes rapidly enough as can be used for the envelope, which we denote below by . By calculating the spectrum of the radiated light from the Fourier transform of the current as a function of time, we can examine the high harmonic generation in the system.

Before jumping into numerics, we first discuss some general principles and examine symmetry considerations. If the lattice is inversion symmetric, then the parity operation immediately tells us that the energy is an even function of the wavevector

 ϵ(k)=ϵ(−k). (8)

If the lattice lacks inversion symmetry, this result still holds if there is no spin-orbit coupling and no magnetic field due to magnetic order. This is because the Hamiltonian is time-reversal symmetric, so there is a Kramers doublet,

 ϵ↑(k)=ϵ↓(−k). (9)

Because the system is also degenerate with respect to spin [since the operator commutes with ], we find , which then implies Eq. (8). Finally, we also examine what happens in the presence of spin-orbit coupling of atomic origin. In this case, since the spin-orbit potential has the periodicity of the lattice, the eigenstates remain in the Bloch form, but since each component of spin no longer commutes with the Hamiltonian, the spin becomes entangled with the wavevector. We continue to have a Kramers doublet, where the spin label in Eq. (9) is replaced by a pseudospin. Because and under a time-reversal operation, the two degenerate states will have opposite velocities. Hence, in equilibrium (), the total current vanishes because the contributions from the states with opposite momentum and opposite pseudospin cancel in the summation in Eq. (5).

Now consider the case where Eq. (8) holds (inversion symmetry or time-reversal invariance and no spin terms in the Hamiltonian). Then, we can replace the exponential in Eq. (1) by a cosine, so we have

 ϵ(k)=−∑δtδcos(k⋅δ). (10)

The field-dependent velocity then becomes

 vk(t)=−1ℏ∑δδtδsin([k−eA(t)ℏc]⋅δ). (11)

Using the trigonometric summation formula, we can expand the sine function into the difference of products of sines and cosines with arguments of and . Substitute this into the equation for the current [Eq. (5)] and recall that the bandstructure is an even function of . This implies that the only term that survives in the summation is the term proportional to , since the sine term is an odd function of the wavevector. This then yields

 J(t)=eℏ∑k∑δδtδcos(k⋅δ)sin(eA(t)⋅δℏc)f(ϵk−μ). (12)

After performing the summation over momentum, the current will be expressed as a constant amplitude for each , , multiplied by the sine function which contains all of the time dependence, and then summed over the neighbors. We manipulate the sine function as follows:

 sin(eA(t)⋅δℏc) = −sin[E0⋅δcos(ω0t+ϕ)g(t)] (13) = −ImeiE0⋅δcos(ω0t+ϕ)g(t) (14) = −Im{J0[E0⋅δg(t)] (16) + 2∞∑n=1inJn[E0⋅δg(t)]cos[n(ω0t+ϕ)]} = 2∞∑n=1(−1)nJ2n+1[E0⋅δg(t)] ×cos[(2n+1)(ω0t+ϕ)].

Here, we use the general formula for the slowly varying envelope function, and the symbol denotes the th Bessel function. The Fourier transform of this will correspond to sharp peaks around the odd harmonics , with the width of the peak determined by how rapidly the envelope function decays in time. The amplitude of the peak can decrease rather rapidly, as the Bessel functions decay for large index, with a fixed argument. Hence, one expects to generically see the odd harmonics in the high harmonic generation if the system obeys the relationship in Eq. (8). It is straightforward to generalize this argument to non-inversion-symmetric systems with spin-orbit coupling if the total current is simply the sum over contributions from the two bands with different pseudospin.

## Iii Numerical results on model systems

We will examine the generic current response for hypercubic lattices in -dimensions, where the field has equal components along each coordinate direction, or has no component along a coordinate direction [hence we examine (1,0,0), (1,1,0), or (1,1,1) directions for a simple cubic lattice in ]. Using the fact that all coordinate directions are equivalent, so is independent of the coordinate direction (with the lattice constant), then gives that the components of the current satisfy

 Ji(t)=ea2ℏd∑δsinAi(t)∫dϵϵρ(ϵ)f(ϵ−μ), (17)

where is the noninteracting density of states of the -dimensional hypercubic lattice. This form of the current involves a temperature-dependent amplitude term and the time-dependent term

 Ji(t)B=sin(E0iℏceacos(ω0t+ϕ)1√πe−(t/t0)2). (18)

Here we define which is the dimensionless amplitude of the electric field (the units of the field are ). The radiated light is then proportional to , which we find by taking the Fourier transform of Eq. (18). Note that this specific form of the radiated harmonics is independent of both the direction of the field and the dimensionality of the lattice for noninteracting electrons on hypercubic lattices.

We begin by showing the normalized electric field in the top panel of Fig. 1. We choose , , and . This produces an even (in time) field profile that has on the order of 10-12 oscillations clearly visible, allowing for the frequency of the pump pulse to be reasonably well defined. The current as a function of time (without the temperature and dimension-dependent amplitude factor ) is shown in the bottom panel. One can clearly see that as the field amplitude rises, the periods of the oscillations in the current get shorter, and that there is a “clipping” like effect in the oscillations when the field amplitude is large enough because the sine function is bounded between plus and minus one. The current has significant structure in it as a function of time and hence will give rise to a complicated Fourier transform. But, by just looking at the results in Fig. 1, we can reason that we expect to see higher frequencies in the Fourier transform when the field amplitude is higher. Indeed, this is what we next show.

In Fig. 2, we plot the logarithm of the radiated light spectrum in arbitrary units, as a function of the frequency of the light. The field profile is identical to that found in Fig. 1, and the different panels show the effect of increasing the amplitude of the pump pulse. The plot is a semilogarithmic plot, because the higher harmonics are often suppressed by orders of magnitude from the fundamental before we hit a “noise floor” in the data. Even though we have a completely well defined integrand to integrate over, and we use an adaptive integration routine available from quadpack for handling such Fourier integrals, it still is challenging to obtain high accuracy for these integrations. The main source of error comes from balancing the integration range versus the order of the polynomials used in the subdivided integration ranges and the absolute accuracy requested in the subroutine. We believe that most of our results with an amplitude larger than are completely accurate, but we do not know whether there might be additional structure in the noise region. Note that this is purely a numerical issue, and will have little bearing on experiment, because once the higher harmonics are too small in amplitude, they cannot be measured. In general, we see on the order of a few up to maybe 20 higher harmonics appearing clearly in the figures.

One can see two main trends occurring as the amplitude is increased: (i) first, the start of the noise floor moves to higher frequencies as the field amplitude is increased, as we expected, and hence we can see more higher harmonics as the system is driven harder, and (ii) we see first the development of a flat plateau of high amplitude harmonics occur at lower frequencies, which increases in size as the field amplitude increases. But, we also see the development of additional structure that looks like noise for these high field-amplitude curves. To investigate this further, we have enlarged the low-frequency region of panel 2 (d) in Fig. 3 for the case with . One can see that this fine structure is not noise, as it appears in the figures with the wider frequency range, but is instead a fine structure of additional oscillations in the harmonic peaks that show more structure the lower the harmonic frequency is. The structure is reminiscent of a beating effect, but we have not been able to completely track down its origin.

In Fig. 4, we investigate the effect of the phase of the field on the high harmonic generation. The effect is rather mild, as can be seen in the figure, where the different colors correspond to different phases, but the evolution is definitely not monotonic in the phase. Since detailed field profiles and phases of the electric field are difficult to determine for a given pump pulse, this weak phase dependence makes the experimental results much more robust, since the most important aspects of the field are the pulse width, oscillation frequency, and amplitude.

## Iv Numerical results for Zinc Oxide

We now apply this model to ZnO using the band structure calculated within density functional theory (DFT). The wurtzite crystal structure of ZnO is uniaxial and non-centrosymmetric, and spin-orbit coupling only weakly affects the band structure (the splitting between bands at the top of the valence band is only a few meV). Thus based on the symmetry arguments presented in Sec. II, we expect to see only odd harmonics in the current response.

The DFT calculations were performed using the Vienna Ab Initio Simulation Package (VASP)  vasp (), with the electron-ion interaction described using projector augmented wave potentials PAW (), and the exchange-correlation interaction handled with the generalized gradient approximation (GGA) PBE (). Spin-orbit effects were not included. The structural parameters were optimized, and the Kohn-Sham eigenvalues were tabulated on uniform grids in the Brillouin zone. Band velocities were calculated using centered differences. We have found that the noise floor in the calculated spectrum converges slowly with respect to the Brillouin-zone sampling. With a grid of k-points, the noise floor in lies about 7 orders of magnitude below the maximum, which is sufficient to distinguish on the order of ten or more harmonics, depending on the magnitude of the field.

The calculated band structure is shown in Fig.  5. The highest lying valence bands have mostly O 2 character, while the conduction band is primarily of Zn 4 character. The GGA severely underestimates the ZnO band gap, yielding a value of about 0.8 eV as compared to the measured gap of 3.4 eV. Quasiparticle calculations based on the GW approximation have found that the main correction to the DFT bands for ZnO is an upward shift by about 2 eV of the two lowest conduction bands vanSchilfgaarde (). The overall shape and width of the O 2 valence bands and of the Zn 4 conduction bands do not change significantly. Since we address here only the contribution to high-order harmonic generation from field-induced acceleration of electrons within bands, and ignore contributions from field-induced dipole transitions between bands, the GGA bands provide a reasonable starting point for a band-by-band analysis.

Because we are not modeling the photoexcitation process, we need to artificially shift the chemical potential to create carriers in this semiconducting system. In the results presented below, the response of electrons in the conduction band is calculated by setting the chemical potential to the equilibrium value corresponding to a 10% occupation of the lowest conduction band. Given the large width of the Zn 4 band, this corresponds to a chemical potential about 2 eV above the conduction band minimum, as marked in Fig.  5. The O 2 bands have larger multiplicity and are much narrower than the conduction band, so a 10% hole occupancy of valence bands requires a much smaller shift in the chemical potential with respect to the band edge. Without a mechanism for interband transitions, the response from holes is estimated by simply summing the current from all the valence bands with partial occupation.

A pulse of the form given by Eq. (8) is applied, with phase , width = 40 fs, frequency , and maximum electric field strengths on the order of 0.1 to 1 V/Å. These parameters yield pump pulses similar to the cycle, mid-infrared laser pulses used in the recent experimental study of high-harmonic generation in ZnO zno_exp (). Figure  6 shows the spectrum of light emitted by electrons in the conduction band when driven by an electric field polarized along the direction. The maximum field strength of the pulse doubles between successive panels, showing the effect of field amplitude on the spectrum. The same general trends that were observed in the calculation based on tight-binding bands are evident here. In particular, when the the field amplitude is increased, a plateau develops in the low-frequency regime, and the low-order harmonic peaks acquire more and more fine structure. In the experimental setup in Ref. zno_exp () the detection cut-off at low-energies was set approximately to the band gap, so the reported spectra start from about the 9th order harmonic and range up to about the 25th harmonic. Within this range, the calculated decay of the intensity with increasing harmonic order agrees roughly with the measured spectra. It would be interesting to try to extend the measurements to lower-order harmonics to see if the fine structure that appears in our calculations is observed.

Given the hexagonal symmetry of the wurtzite crystal structure, it is expected that the response of electrons in ZnO will depend on whether the electric field is parallel or perpendicular to the axis. In addition, since the dispersion of the conduction and valence bands are quite different, it is interesting to compare the emission spectra for electrons and holes. Figure  7 shows the spectra generated by electrons in the conduction band and by holes in the valence band when driven by fields in the and directions. The parameters for the pump pulse (with the exception of the field direction) are the same as those used for Fig.  6c. When the driving field lies in the basal plane (Fig.  7b and d), the low-frequency plateau region is shortened compared to when the field is perpendicular to the plane (Fig. 7a and c). The in-plane field produces a stronger signal for the lower harmonics, but the magnitude of the harmonic peaks decays more quickly with harmonic order. Further, the low-order peaks are narrower and exhibit less structure than when the field is polarized perpendicular to the plane. Despite significant differences in for the electrons and holes, the spectral characteristics of their radiation are very similar. There is an overall reduction in signal intensity from the valence band as compared to the conduction band, which is consistent with the holes having a signficantly larger band mass than the electrons. Interestingly, the fine structure in the lower harmonics seems to depend much more on the direction of the electric field than on whether the carriers are in the valence or conduction band. At higher frequencies, the detailed shape of the peaks exhibits stronger dependence on the band structure. Varying the filling factors for the conduction and valence bands changes some of the detailed features in the spectra as well as the overall intensity of emitted light, but the general trends discussed here seem to be generic and robust.

A remaining puzzle is that while most crystallographic orientations produced only odd-order harmonics in the experiments, both even and odd harmonics were observed when the polarization of the incident electric field was parallel (or close to parallel) to the axis of the ZnO crystal zno_exp (). Our model based on intraband dynamics of carriers in the bulk is not able to account for the appearance of even harmonics under these circumstances. To more fully understand the experimental results for ZnO (or other semiconductors), the coupling between bands golde () and other effects such as the broken symmetry at the surface need to be further explored.

## V Conclusions

In this work, we have shown some simple theories for high harmonic generation in metals (or photodoped semiconductors) which capture many of the experimental results already seen in this process. While it has long been thought that Bloch oscillations occur at too fast a time scale to be observed experimentally in metals, this particular experiment is the exception, which allows one to indirectly observe the oscillations from the light generated by the accelerating electrons. Of course, in any real system, one needs to include interactions to determine how feasible the observation of these effects are. The most important ones will be interactions with defects (disorder), with phonons, and with other electrons as the correlations in the material grow. In order to do this, a more complete many-body physics formulation is needed on the Kadanoff-Baym-Keldysh contour kadanoff_baym (); keldysh (), which we have carried out an initial study for elsewhere kemper (). It will be interesting to examine how this behavior survives as correlations increase, into materials like Mott insulators. In addition, it will be important to understand the balance between the more universal behavior, which is captured by model system calculations, and the materials specific features, which require more materials specific calculations based on density functional theory plus extensions to include correlations in the nonequilibrium domain. The future for these types of studies, both theoretical and experimental is clearly a bright one!

###### Acknowledgements.
The DFT calculations (AYL and JKF) were supported by the National Science Foundation under grant number DMR-1006605. The model system calculations were supported by the U.S. Department of Energy, Basic Energy Sciences, Materials Sciences and Engineering Division under contract No. DE-FG02-08ER46542 (JKF) and contract No. DE-AC02-76SF00515 (AFK and TPD). JKF was also supported by the McDevitt bequest at Georgetown University. The Georgetown-Stanford collaboration was supported by the U.S. Department of Energy, Basic Energy Sciences, Materials Sciences and Engineering Division under contract No. DE-FG02-08ER46540.

### References

1. F. Bloch, Z. Phys. 52, 555 (1928).
2. C. Zener, Proc. R. Soc. A 145, 523 (1934).
3. C. Waschke, H. G. Roskos, R. Schwedler, K.l Leo, and H. Kurz Phys. Rev. Lett. 70, 3319 (1993); V. Lyssenko, G. Valuŝis, F. Löser, and T. Hasche, Phys. Rev. Lett. 79, 301 (1997).
4. B. M. Dahan, E. Peik, J. Reichel, Y. Castin, and C. Salomon, Phys. Rev. Lett. 76, 4508 (1996).
5. A. F. Kemper, B. Moritz, J. K. Freericks, and T. P. Devereaux, unpublished, arxiv:1204.1803.
6. D. Golde, T. Meier, and S. W. Koch, Phys. Rev. B 77, 075330 (2008).
7. S. Ghimire, A. D. DiChiara, E. Sistrunk, P. Agostini, L. F. DiMauro, and David A. Reis, Nat. Phys. 7, 138 (2011).
8. V. M. Turkowski and J. K. Freericks, Phys. Rev. B 71, 085104 (2005).
9. L. P. Kadanoff and G. Baym, Quantum Statistical Mechanics, Benjamin, New York, 1962.
10. L. V. Keldysh, Zh. Eksp. Teor. Fiz. 47, 1515 (1964).
11. G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996); Phys. Rev. B 54, 11169 (1996).
12. P. E. Blochl, Phys. Rev. B 50, 17953 (1994); G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
13. J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
14. M. Usuda, N. Hamada, T. Kotani, and M. van Schilfgaarde, Phys. Rev. B 66, 125101 (2002).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters