Gauge invariance and the interpretation of inter- and intraband processes in high-order harmonic generation from bulk solids

TL;DR

This paper presents a gauge-invariant theoretical model for high-order harmonic generation in bulk solids, emphasizing the equal treatment of inter- and intraband currents and highlighting the gauge dependence of their interpretation.

Contribution

It introduces a gauge-invariant approach to modeling HHG in solids, clarifying the roles of inter- and intraband processes and their interpretation.

Findings

Total HHG source is gauge independent

Inter- and intraband distinctions are gauge dependent

Model aligns with both analytic and numerical results

Abstract

A theoretical model for high-order harmonic generation (HHG) in bulk solids is considered. Our approach treats laser-induced inter- and intraband currents on an equal footing. The sum of these currents is the source of the high-order harmonic radiation, and does not depend on the particular electromagnetic gauge we choose to describe the process. On the other hand, as it is shown using analytic and numerical calculations, the distinction between intra- and interband dynamics is gauge dependent, implying that the interpretation of the process of HHG using these terms requires carefulness.

Equations

$$\hslash \frac{d\mathbf{k}(t)}{dt}=e_0\mathbf{F}(t),$$

$$H(t)=$$

$$=$$

$$H_{ext}^v(t)=\frac{1}{2m}\left[-2e_0\mathbf{p}{\mathbf{A}}^v(t)+e_0^2{A}^{v2}(t)\right].$$

$$H_{ext}^l(t)=-e_0\mathbf{F}(t)\mathbf{r}.$$

$${\Psi }^l(\mathbf{r},\mathbf{t})=U(\mathbf{r},t){\Psi }^v(\mathbf{r},t)=\exp \left[-\frac{i}{\hslash }{e}_0\mathbf{A}(t)\mathbf{r}\right]{\Psi }^v(\mathbf{r},t).$$

$$|{\Psi }^v(t)⟩=\sum _{n\mathbf{k}}{c}_{n\mathbf{k}}(t)|n,\mathbf{k}⟩,|{\Psi }^l(t)⟩=\sum _{n\mathbf{k}}{b}_{n\mathbf{k}}(t)|n,\mathbf{k}⟩.$$

$$b_{n^{\prime }{\mathbf{k}}^{\prime }}(t)=\sum _{n\mathbf{k}}⟨n^{\prime }{\mathbf{k}}^{\prime }|\exp \left[-i{e}_0\mathbf{A}(t)r/\hslash \right]|n,\mathbf{k}⟩c_{n\mathbf{k}}(t).$$

$$\rho (t=0)=\sum _{n\mathbf{k}}|n,\mathbf{k}⟩⟨n,\mathbf{k}|f[E_n(\mathbf{k})],$$

$${\rho }^v(t)=\sum _{n{n}^{\prime }\mathbf{k}}{c}_{n\mathbf{k}}(t)c_{n^{\prime }\mathbf{k}}^{\ast }(t)|n,\mathbf{k}⟩⟨n^{\prime },\mathbf{k}|,$$

$$\mathbf{J}:=\frac{e_0}{\mathcal{V}m}\mathrm{Tr}\left[{\rho }^g{\mathbf{p}}_{kin}\right]=\frac{e_0}{\mathcal{V}m}\mathrm{Tr}\left[{\rho }^g(\mathbf{p}-e_0{\mathbf{A}}^g)\right].$$

$$\mathbf{J}$$

$$-i\hslash \frac{e_0}{\mathcal{V}m}\sum _{n,n^{\prime }}\sum _{\mathbf{k}}{c}_{n\mathbf{k}}^{\ast }(t)c_{n^{\prime }\mathbf{k}}(t){\int }_{\Omega }{u}_{n,\mathbf{k}}^{\ast }(\mathbf{r})\nabla u_{n^{\prime },\mathbf{k}}(\mathbf{r})d^3\mathbf{r}.$$

$$\mathbf{J}$$

$$-i\hslash \frac{e_0}{\mathcal{V}m}\sum _{n,n^{\prime },\mathbf{k}}{b}_{n\mathbf{k}}^{\ast }(t)b_{n^{\prime }\mathbf{k}}(t){\int }_{\Omega }{u}_{n,\mathbf{k}}^{\ast }(\mathbf{r})\nabla u_{n^{\prime },\mathbf{k}}(\mathbf{r})d^3\mathbf{r}.$$

$${\mathbf{j}}^g=\frac{e_0}{\mathcal{V}m}\sum _{n\mathbf{k}}{\left|c_{n\mathbf{k}}^g(t)\right|}^2⟨n,\mathbf{k}|{\mathbf{p}}_{kin}|n,\mathbf{k}⟩,$$

$${\dot{\mathbf{P}}}^g=\frac{e_0}{\mathcal{V}m}\sum _{n\ne n^{\prime }\mathbf{k}}{\left[c_{n\mathbf{k}}^g(t)\right]}^{\ast }{c}_{n^{\prime }\mathbf{k}}^g(t)⟨n,\mathbf{k}|{\mathbf{p}}_{kin}|n^{\prime },\mathbf{k}⟩$$

$$b_{n^{\prime }{\mathbf{k}}^{\prime }}(t)=\sum _{n\mathbf{k}}⟨n^{\prime }{\mathbf{k}}^{\prime }|\exp \left[-i{e}_0\mathbf{A}(t)r/\hslash \right]|n,\mathbf{k}⟩c_{n\mathbf{k}}(t)$$

$$=\sum _{n\mathbf{k}}{\int }_{\mathcal{V}}{e}^{i\mathbf{r}(\mathbf{k}-e\mathbf{A}(t)/\hslash +{\mathbf{k}}^{\prime })}{u}_{n^{\prime },{\mathbf{k}}^{\prime }}^{\ast }(\mathbf{r})u_{n,\mathbf{k}}(\mathbf{r})d^3\mathbf{r}{c}_{n\mathbf{k}}(t)$$

$$=\sum _n{\int }_{\Omega }{u}_{n^{\prime },{\mathbf{k}}^{\prime }}^{\ast }(\mathbf{r})u_{n,\mathbf{k}(t)}(\mathbf{r})d^3\mathbf{r}{c}_{n\mathbf{k}(t)}(t),$$

$${\mathbf{j}}^v=\frac{e_0}{\mathcal{V}m}\sum _{n\mathbf{k}}{\left|c_{n\mathbf{k}}(t)\right|}^2⟨n,\mathbf{k}|\mathbf{p}|n,\mathbf{k}⟩-\frac{\mathcal{N}{e}_0^2}{\mathcal{V}m}\mathbf{A},$$

$${\mathbf{j}}^l=\frac{e_0}{\mathcal{V}m}\sum _{n\mathbf{k}}{\left|b_{n\mathbf{k}}(t)\right|}^2⟨n,\mathbf{k}|\mathbf{p}|n,\mathbf{k}⟩.$$

$$⟨n,\mathbf{k}|\mathbf{p}|n^{\prime },{\mathbf{k}}^{\prime }⟩=$$

$$=-i\hslash {\int }_{\mathcal{V}}{\Psi }_{n,\mathbf{k}}^{\ast }(\mathbf{r})\nabla {\Psi }_{n^{\prime },{\mathbf{k}}^{\prime }}(\mathbf{r})d^3\mathbf{r}$$

$$=-\frac{i\hslash }{\mathcal{V}}{\int }_{\mathcal{V}}{e}^{i\mathbf{r}({\mathbf{k}}^{\prime }-\mathbf{k})}{u}_{n,\mathbf{k}}^{\ast }(\mathbf{r})(i{\mathbf{k}}^{\prime }+\nabla )u_{n^{\prime },{\mathbf{k}}^{\prime }}(\mathbf{r})d^3\mathbf{r}$$

$$=-\frac{i\hslash }{\mathcal{V}}\sum _{\mathbf{R}}{e}^{i\mathbf{R}({\mathbf{k}}^{\prime }-\mathbf{k})}{\int }_{\Omega }{u}_{n,\mathbf{k}}^{\ast }(\mathbf{r})(i{\mathbf{k}}^{\prime }+\nabla )u_{n^{\prime },{\mathbf{k}}^{\prime }}(\mathbf{r})d^3\mathbf{r}$$

$$={\delta }_{\mathbf{k}{\mathbf{k}}^{\prime }}\left[\hslash \mathbf{k}{\delta }_{n{n}^{\prime }}-i\hslash {\int }_{\Omega }{u}_{n,\mathbf{k}}^{\ast }(\mathbf{r})\nabla u_{n^{\prime },{\mathbf{k}}^{\prime }}(\mathbf{r})d^3\mathbf{r}\right],$$

$$⟨n,\mathbf{k}|\frac{1}{2m}\left[-2e_0\mathbf{pA}(t)+e_0^2{A}^2(t)\right]|n^{\prime },{\mathbf{k}}^{\prime }⟩=$$

$$=-\frac{e_0}{m}⟨n,\mathbf{k}|\mathbf{p}|n^{\prime },{\mathbf{k}}^{\prime }⟩\mathbf{A}(t)+\frac{1}{2m}{e}_0^2{A}^2(t){\delta }_{\mathbf{k}{\mathbf{k}}^{\prime }}{\delta }_{n{n}^{\prime }}$$

$$={\delta }_{\mathbf{k}{\mathbf{k}}^{\prime }}\left[-i\hslash {\int }_{\Omega }{u}_{n,\mathbf{k}}^{\ast }(\mathbf{r})\nabla u_{n^{\prime },{\mathbf{k}}^{\prime }}(\mathbf{r})d^3\mathbf{r}\right]\mathbf{A}(t)$$

$$+\left[\hslash \mathbf{k}\mathbf{A}(t)+\frac{1}{2m}{e}_0^2{A}^2(t)\right]{\delta }_{\mathbf{k}{\mathbf{k}}^{\prime }}{\delta }_{n{n}^{\prime }}.$$

$$\left[\frac{{\hslash }^2}{2m}{\left(-i\frac{\partial }{\partial x}+k\right)}^2+V(x)\right]u_{n,\mathbf{k}}(\mathbf{r})=E_n(k)u_{n,\mathbf{k}}(\mathbf{r}).$$

$$S_{n{n}^{\prime }}^{k{k}^{\prime }}={\int }_{-a/2}^{a/2}{u}_{n,k}^{\ast }(x)u_{n^{\prime },k^{\prime }}(x)dx$$

$$A(t)=A_0{\sin }^2\left(\frac{\pi t}{T}\right)\cos (\omega t),$$

Full text

*Introduction * The idea that solid state targets illuminated by strong, short near infrared pulses can produce coherent X-ray radiation Kálmán and Brabec (1995) or even attosecond electromagnetic bursts Faisal and Kamiński (1996) preceded the first experiments that demonstrated the appearance of high-order harmonics (up to order of 25) in the spectra of strongly driven solid state samples. In Ref. Ghimire et al. (2011) a wide bandgap ZnO target (3.2 eV) was excited by a pulse with central wavelength of 3.25 $\mu$m, meaning that at least 9 photons are required for an excitation from the valance to the conduction band [see also Ghimire et al. (2014) for more details]. Recently, a direct comparison of high-order harmonic generation in the solid and gas phases of argon and krypton has been reported Ndabashimiye et al. (2016).

Theoretically, semiconductor Bloch-equations were applied to describe the problem Golde et al. (2008), a closed-form expression were given to the subcycle-resolved transition rate of electrons between bands Hawkins and Ivanov (2013), appearance of attosecond pulses were predicted Higuchi et al. (2014), semiclassical Vampa et al. (2014) and a saddle-point Vampa et al. (2015) analysis were performed and the role of an indirect bandgap was also investigated Yamaguchi and Nobusada (2016).

The intensity of the exciting laser pulse was below the damage threshold in Ref. Ghimire et al. (2011). Although high-order harmonic generation (HHG) is known to be possible also in plasmas emerging from solid state surfaces as a consequence of intense electromagnetic radiation [see e.g., Teubner and Gibbon (2009) for a review], here we focus on the case when a system of electrons in a static, periodic potential interacts with the exciting laser field Hawkins et al. (2015).

Neglecting relaxation and multiparticle effects, the response of the electronic subsystem in a solid to the external laser pulse is usually categorized as follows: (1) Interband transitions, (2) Laser-driven intraband motion, (3) Dynamical Bloch-oscillations. The first phenomenon is the usual optical generation of charge carriers, but this is already a nonlinear effect in the range of HHG. Points 2 and 3 above are closely related, it is only the strength of the laser field that determines whether dynamical Bloch oscillations appear or not. As a quasiclassical picture, the external field of the laser changes the crystal momentum of the electrons. More precisely, as described by the acceleration theorem Kittel (2005); Arlinghaus and Holthaus (2011), the change of the crystal momentum can be written as

[TABLE]

where $\mathbf{F}(t)$ represents the time dependent external electric field, and $e_{0}$ denotes the charge of the electron. (That is, here, and in the following, $e_{0}$ is negative.) When the amplitude of this oscillatory motion is large enough to cross the boundary of the first Brillouin zone, Bragg-reflection occurs, which can be termed as dynamical Bloch-oscillation Mücke (2011); Földi et al. (2013). It is important to note that mechanisms 1-3 are coupled, practically none of them appears on its own Higuchi et al. (2014).

As we show in the following, although separating inter- and intraband dynamics is very useful for the intuitive interpretation of the physical processes, it depends on the choice of the electromagnetic gauge. In other words, gauge transformations (which, obviously, do not change physically relevant results) mix the inter- and intraband dynamics.

The Hamiltonian describing the electron in a periodic potential (representing the solid) and interacting with an external field can be written as:

[TABLE]

where $V(\mathbf{r})$ is the lattice-periodic potential of the crystal, while $H_{ext}(t)$ takes the interaction with the external field into account, via the electromagnetic potentials. Note the appearance of the kinetic momentum $\mathbf{p}_{kin}=m\mathbf{v=p-}e_{0}\mathbf{A}$ in Eq. (2), which is seen to be different from the canonical $\mathbf{p}=-i\hbar\nabla$  depending on the gauge.

The eigenstates of the field-free Hamiltonian obeying $H_{0}|n,\mathbf{k}\rangle=E_{n}(\mathbf{k})|n,\mathbf{k}\rangle,$ in coordinate representation, take the form $\Psi_{n,\mathbf{k}}(\mathbf{r})=\exp(i\mathbf{kr})u_{n,\mathbf{k}}(\mathbf{r})/\sqrt{\mathcal{V}},$ where $n$ denotes the band index of these Bloch-states, $u_{n,\mathbf{k}}(\mathbf{r})$ are lattice periodic functions and $\mathcal{V}$ is the crystal volume. We assume here, as usual, periodic boundary conditions, and then the $\mathbf{k}$-space will not be continuous, we obtain a discrete, densely spaced series of $\mathbf{k}$ vectors. The $\mathbf{k}$-dependent eigenenergies $E_{n}(\mathbf{k})$ provide the dispersion relations for the various bands.

In order to specify the interaction term $H_{ext}(t)$, we have to choose a gauge. Clearly, physical predictions based on exact solutions without approximations must be the same in all gauges. We assume here that the spatial variation of the external field $\mathbf{F(}t\mathbf{)}$ can be neglected, since the wavelengths of the exciting infrared pulses are much longer than the lattice constants in a crystal. Then one of the appropriate gauges is the velocity ($v$) gauge, where the vector and scalar potentials of the field $\mathbf{F}(t)$ are $\mathbf{A}^{v}(t)=-\int_{-\infty}^{t}\mathbf{F}(t^{\prime})dt^{\prime},\Phi^{v}=0$. With this choice, we have

[TABLE]

The second widely used option is the length ($l$) gauge, in which $\mathbf{A}^{l}=0,\Phi^{l}=-\mathbf{rF}(t)$, obtained from the $v$ gauge by using the gauge function $\Lambda=-\mathbf{rA}^{v}(t)$. This yields

[TABLE]

(Note that both of these choices belong to the class of Coulomb gauges.) From now on, we shall use the notation $\mathbf{A}(t)\equiv\mathbf{A}^{v}(t)$. It should be kept in mind that electron wavefunctions in these gauges are connected by a space-time dependent unitary transfomation:

[TABLE]

In the following we will consider pulsed excitation of duration $T,$ and assume that both $\mathbf{F}(t)$ and $\mathbf{A}(t)$ vanish outside the interval of $[0,T].$ This means that states in the two gauges coincide for $t<0$ and $t>T$ [when $U(\mathbf{r},t)$ is the identity]. Note that this choice avoids gauge-related ambiguities for $t>T$ Madsen (2002); Rzazewski and Boyd (2004); Han and Madsen (2010).

Under the effect of $H_{ext}^{g},$ any initial state will evolve in time, and can be expanded as a time dependent linear combination of the unperturbed eigenstates: As $H_{ext}^{g}$ is gauge dependent, the expansion coefficients as well as the states $\left|\Psi^{g}(t)\right\rangle$ shall depend on the gauge chosen. We use the following notations in the two gauges of interest here:

[TABLE]

The following (invertible) relation holds between the expansion coefficients:

[TABLE]

Before the arrival of the exciting laser pulse ($t<0$), the solid state target can be assumed to be in thermal equilibrium. This initial condition cannot be described by a pure quantum mechanical state, we have to consider a (single particle) density operator, which is diagonal in the eigenstates of the unperturbed Hamiltonian $H_{0}:$

[TABLE]

where the relative statistical weights of the projectors is determined by the Fermi function $f$. For wide-bandgap target materials at room temperature, $f[E_{n}(\mathbf{k})]$ is practically unity for the valence band, and zero for the conduction bands. That is, to a very good approximation, we can write $\rho(t=0)=\sum_{\mathbf{k}}|n_{0},\mathbf{k}\rangle\langle n_{0},\mathbf{k}|,$ where $n_{0}$ corresponds to the valence band. Similarly to the pure quantum mechanical states, the density operator is also gauge dependent, the analogue of Eq. (5) reads $\rho^{l}(t)=U(t)\rho^{v}(t)U^{\dagger}(t)$ and the time evolution of the density operator is given by the von-Neumann equation $i\hbar\frac{\partial}{\partial t}\rho^{g}(t)=\left[H^{g}(t),\rho^{g}(t)\right].$ For the sake of simplicity, let us use the velocity gauge, when (in dipole approximation) the $A^{2}$ term in the Hamiltonian commutes with $\rho^{v}$, and the remaining part does not mix states with different indices $\mathbf{k}$ (see the Appendix for more details). Therefore we have:

[TABLE]

where $\sum_{n}|c_{n\mathbf{k}}(t)|^{2}=1$, and the time dependence of these coefficients is determined by the von-Neumann equation [with the initial conditions of $c_{n\mathbf{k}}(0)=\delta_{n,n_{0}}$]. Since the trace of the density operator is preserved during the dynamics, our choice of normalization means that $\mathrm{Tr}\rho$ equals the number $\mathcal{N}$ of valence band states in the first Brillouin zone. Note that this approach does not provide absolute field strengths, exact ”number of photons” in the HHG modes (unlike Ref. Gombkötő et al. (2016)), it is only the relative weight of the high harmonics that can be obtained. That is the reason why we are free to choose the normalization condition for the initial density matrix.

*High-order harmonic generation * Moving charges mean the source of the high harmonic radiation, in other words it is the expectation value of the time dependent current (density) that is to be calculated:

[TABLE]

It is important to stress here, that the source of any secondary radiation (such as high harmonics) is not the canonical, but the kinetic momentum $\mathbf{p}_{kin},$ due to its direct connection to the velocity operator Cohen-Tannoudji et al. (2004).

Let us look at the expectation value above in the two gauges we consider here. Similarly to the more detailed calculations presented in the Appendix, in the velocity gauge we obtain

[TABLE]

In the length gauge $\mathbf{A}^{l}\mathbf{=}0$, and we have therefore:

[TABLE]

The results calculated in any of the gauges must be identical, but this equality is valid only for the whole expressions in Eqs. (11) and (12), and *by no means for the individual terms on the right hand sides. *

*Gauge dependence * It is tempting to regroup the terms in Eqs. (11) and (12) and interpret the macroscopic, ”current-like” part of the source of the HHG radiation

[TABLE]

as the intraband contribution, while the microscopic, ”polarization-like” part

[TABLE]

as the interband contribution. Indeed, $\langle n,\mathbf{k}|\mathbf{p}_{kin}|n,\mathbf{k}\rangle/m$ is the band velocity in the $n$th band at point $\mathbf{k},$ thus $\mathbf{j}^{(g)}$ means the sum of all these velocities weighted by the corresponding populations. For comparison with textbook methods Haug and Koch (2004), let us note that $\langle n,\mathbf{k}|\mathbf{p}|n,\mathbf{k}\rangle$ can also be calculated as $m/\hbar\nabla_{\mathbf{k}}E_{n}(\mathbf{k}),$ but Eqs. (11) and (12) treat inter- and intraband currents on an equal footing. Clearly, $\mathbf{J}=\mathbf{j}^{g}+\dot{\mathbf{P}}^{g},$ and it is gauge independent, but it does not necessarily hold for $\mathbf{j}^{g}$ and $\dot{\mathbf{P}}^{g}$ separately.

Let us investigate the relation (7) between the expansion coefficients $c_{n\mathbf{k}}(t)$ and $b_{n\mathbf{k}}(t).$ Assume that the dynamics is solved in the velocity gauge [i.e., $c_{n\mathbf{k}}(t)$ are known functions of time], and let us transform the result to length gauge. We obtain:

[TABLE]

where $\mathbf{k}(t)=e\mathbf{A}(t)/\hbar-\mathbf{k}^{\prime}.$ The first consequence of this result is that whenever the velocity gauge density operator is diagonal in $\mathbf{k},$ [which is the case during the whole time evolution for the initial conditions given by Eq. (8)], the same holds also for the length gauge, but with $\mathbf{k}$ values being shifted by $e\mathbf{A}(t)/\hbar$. [This is in agreement with the acceleration theorem (1).] On the other hand, since the lattice-periodic functions $u_{n,\mathbf{k}}$ are orthogonal only for the same index $\mathbf{k},$ the integral in the last line of Eq. (15) is not proportional to $\delta_{nn^{\prime}}.$

This effect is closely related to the fact that the separation of $\mathbf{J}$ as the sum of intraband ($\mathbf{j}$) and interband ($\dot{\mathbf{P}}$) components is not unambiguous. E.g., for $\mathbf{j}$ we have

[TABLE]

where $\mathcal{N}=\mathrm{Tr}\rho.$ On the other hand,

[TABLE]

The last term in Eq. (16) is determined solely by the time dependence of the vector potential, while the coefficients $c_{n\mathbf{k}}$ and $b_{n\mathbf{k}}$ that mean the difference between the sums in Eqs. (16) and (17) are related via integrals that are eventually determined by the structure of the material. This means that even if for a certain material and time dependent vector potential $\mathbf{j}^{l}(t)$ and $\mathbf{j}^{v}(t)$ may happen to be identical, there is no physical reason why this should hold e.g., for a different target material. In other words, the inter- and intraband contributions that constitute the gauge independent source term $\mathbf{J},$ cannot be determined without referring to a particular gauge.

*Numerical example *

As an illustration, we calculate the current generated by a short electromagnetic pulse in a one-dimensional model crystal, where the exciting field is polarized along the periodic chain. We use a model potential that produces a direct bandgap of 3.2 eV, similarly to the case of ZnO Ghimire et al. (2011). Having determined the Bloch-sates, the time dependence of $j^{g}$ and $\dot{P}^{g}$ [the one-dimensional versions of the quantities given by Eqs. (13) and (14), respectively] can be calculated, as well as their sum, $J^{g}.$ [The superscript reminds us the fact that various approximations may result in the (incorrect) gauge dependence of $J.$]

As we can see in the bottom panel of Fig. 1, the time dependence of $j^{v}$ and $j^{l}$ are clearly different. The total currents $J^{v}$ and $J^{l}$ are not exactly the same either, but for the 8-band model presented in Fig. 1, the simple one-dimensional numerical model supports our earlier findings within reasonable precision.

Discussion An overview of the results above shows that the main reason why one cannot make a gauge-independent distinction between intra- and interband dynamics is related to the fact that states corresponding to different gauges are not identical. As an analogue, let us recall that a similar statement holds also for a single atom in an electromagnetic field. Even in dipole approximation, depending on the choice of the gauge, the states may ”incorporate some interaction with the field” Bauer et al. (2005). However, for optical frequencies and laser field strengths around 1 GV/m, the phase factor appearing in $\exp(-ie\mathbf{A}\mathbf{r}/\hbar)$ is almost constant over the size of an atom, its change is around $\pi/100,$ thus the effect is weak. In the case of solids, on the other hand, the interaction region is much larger, and the effect cannot be neglected.

Finally, let us emphasize, that all the considerations above are relevant during the excitation. When the exciting pulse is over (i.e., $\mathbf{A}=0$), the separation $\mathbf{J}=\mathbf{j}+\dot{\mathbf{P}}$ is clearly unambiguous, and e.g., $\mathbf{j}$ – as a macroscopic current – can be measured.

Summary We investigated the process of high-order harmonic generation in bulk solids. Working in the single electron picture, light-matter interaction were considered in two widely used electromagnetic gauges. Inter- and intraband contributions to the whole current were treated in a unified way. Both analytic and numerical results show that the interpretation of the process in terms of inter- and intraband currents has to performed with care, since these contributions do depend on the choice of the electromagnetic gauge.

Acknowledgments

The author thanks M. G. Benedict, S. Varró, K. Varjú and A. Czirják for useful discussions. Partial support by the ELI-ALPS project is acknowledged. The ELI-ALPS project (GOP-1.1.1-12/B-2012-000, GINOP-2.3.6-15-2015-00001) is supported by the European Union and co-financed by the European Regional Development Fund. The work was also supported by the Hungarian National Research, Development and Innovation Office under Contract No. 124750.

I Appendix

I.1 Bloch-basis matrix elements

For practical calculations, the matrix elements of the Hamiltonian in the Bloch-state basis are very useful. Remaining in the dipole approximation, where the spatial dependence of the field amplitude can be neglected, we obtain:

[TABLE]

where $\mathbf{R}$ denotes lattice vectors. In view of this, the matrix elements of $H_{ext}$ in velocity gauge [recall that $\mathbf{A(}t\mathbf{)\equiv A}^{v}\mathbf{(}t\mathbf{)}$] are given by:

[TABLE]

Note that the reason why these equations have a relatively simple form is the compatibility of the momentum operator with the translational symmetry of the crystal. On the other hand, the matrix elements of $H_{ext}^{l}(t)$ as given by $-e_{0}\langle n,\mathbf{k}|\mathbf{r}|n^{\prime},\mathbf{k}^{\prime}\rangle\mathbf{F}(t)$ cannot be simplified as above, because the operator $\mathbf{r}$ is not lattice-periodic.

I.2 Numerical details

In order to be able to calculate the matrix elements as discussed above, we need the actual Bloch-states, that is, in a discretized picture, the complex values $u_{nk}(x_{m}),$ where $x_{m}$ is in the unit cell ($-a/2<x_{m}<a/2$). For the sake of definiteness, the lattice constant $a$ is chosen to be 0.5 nm. Once the potential is specified, these eigenfunctions can be obtained using

[TABLE]

As numerical examples, we investigated two model potentials. $V_{1}$ contains two localized attracting centers: $V_{1}(x)[eV]=-25\cos^{2}\left[\pi(x-x_{1})/(15a)\right]$ $-25\cos^{2}\left[\pi(x-x_{2})/(15a)\right]$, where $x_{1}/a=-0.2$ and $x_{2}/a=0.107,$ and the $\cos^{2}$ functions are zero unless their argument is in the interval $[-\pi/2,\pi/2]$. The second model potential is a modification of the one that has been used in Ref. Hawkins et al. (2015): $V_{2}[eV]=−25[1+\tanh(x+x_{0})][1+\tanh(−x+x_{0})],$ where $x_{0}=0.2475a_{0}$ ($a_{0}$ denotes the Bohr radius). There is an important point these potentials have in common: both produce a bandgap of 3.2 eV, which is in agreement with the material properties of ZnO. Clearly, apart from the bandgap, the two potentials correspond to different band schemes, but the qualitative results related to the question of gauge invariance are the same in both cases. For the sake of definiteness, we used $V_{1}$ for obtaining the data presented in Fig. 1.

Having $N$ grid points $x_{n}$ in the 1D unit cell, Eq. (20) with periodic boundary conditions provide a set of eigenfunctions $\left\{u_{n,k}\right\}_{1}^{N}$ that forms a basis for each value of $k.$ (That is, any finite-valued periodic function defined on the grid can be written as a linear combination of these eigenfunctions.) It is convenient to use discrete $k$ indices as well, $k_{m}=-\pi/a+(m-1)\Delta k,\ m=1,\ldots,M$ and $(M-1)\Delta k=2\pi/a.$ In this way, we obtain $N\times M$ functions $u_{n,k},$ and – besides the matrix elements (19) – we can calculate their overlaps

[TABLE]

numerically. Note that – by construction – we have $S_{nn^{\prime}}^{kk}=\delta_{nn^{\prime}}.$ With sufficient number of real-space discretization points, the integrals (19) and (21) converge for the lowest lying energy eigenstates. (According to our experience, it is sufficient for both $N,M$ to have the order of magnitude of a few hundred.)

In velocity gauge, the index $k$ of the Bloch-states does not change during the time evolution, and by inspecting the populations we can determine how many bands play observable role in the time evolution. Using the appropriate (converged) matrix elements, we can follow the time evolution, and transform the result to length gauge at appropriate time instants. More precisely, for a given initial index $k_{0},$ whenever $\mathbf{k}(t)=e\mathbf{A}(t)/\hbar-\mathbf{k}_{0}$ equals one of the k-space grid points $k_{m}$, we can use the overlap matrix to calculate the transformation. Since these time instants $t_{m}$ are the same for any initial index $k_{0},$ the time evolution of the complete density matrix can be calculated and compared in both gauges at the discrete $t_{m}$ values. In this way, discretization in k-space affects only the time resolution of the results: at $t_{m},$ the calculated physical quantities are numerically exact in the sense that the only approximation is real space discretization.

In our actual calculations, the exciting laser field is assumed to be polarized along the $x$ direction, and the time dependence of the only nonzero component of the vector potential is given by

[TABLE]

provided $t\in[0,T],$ and zero otherwise. Similarly to Ref. Ghimire et al. (2011), we consider many-cycle, mid-infrared excitation. That is, $\lambda=3\mu$m and $T=300$ fs in the calculations. The amplitude is chosen to correspond to a peak laser field strength of 1.0 GV/m.

Finally, note that the overlap matrix (21) is clearly unitary, provided all the bands are taken into account. However, when we do not use all the bands (either below or above the highest valence band), the appropriate projection of $S$ may not be exactly unitary, leading to the fact that the norm (trace of the density matrix) is not the same in the two gauges. In other words, besides the usual criteria a numerical method has to satisfy (convergence, stability, etc.), the requirement of gauge invariance sets a new one. Working in Bloch-state basis, this simply means that a large enough number of bands has to be taken into account.