{'fraction_non_alphanumeric': 0.04528104575163399, 'fraction_numerical': 0.009359477124183007, 'mean_word_length': 4.387605633802817, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gauge invariance and the interpretation of inter- and intraband processes in high-order harmonic generation from bulk solids

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.

Introduction The idea that solid state targets illuminated by strong, short near infrared pulses can produce coherent X-ray radiation [1] or even attosecond electromagnetic bursts [2] 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. [3] a wide bandgap ZnO target (3.2 eV) was excited by a pulse with central wavelength of 3.25 µm, meaning that at least 9 photons are required for an excitation from the valance to the conduction band [see also [4] 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 [5].
Theoretically, semiconductor Bloch-equations were applied to describe the problem [6], a closed-form expression were given to the subcycle-resolved transition rate of electrons between bands [7], appearance of attosecond pulses were predicted [8], semiclassical [9] and a saddlepoint [10] analysis were performed and the role of an indirect bandgap was also investigated [11].
The intensity of the exciting laser pulse was below the damage threshold in Ref. [3]. 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., [12] 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 [13].
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 [14,15], the change of the crystal momentum can be written as where 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, Braggreflection occurs, which can be termed as dynamical Bloch-oscillation [16,17]. It is important to note that mechanisms 1-3 are coupled, practically none of them appears on its own [8].
As we show in the following, although separating interand 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: where V (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 p kin = mv = p−e 0 A in Eq. (2), which is seen to be different from the canonical p = −i ∇ depending on the gauge. The eigenstates of the field-free Hamiltonian obeying H 0 |n, k = E n (k)|n, k , in coordinate representation, take the form Ψ n,k (r) = exp(ikr)u n,k (r)/ √ V, where n denotes the band index of these Bloch-states, u n,k (r) are lattice periodic functions and V is the crystal volume. We assume here, as usual, periodic boundary conditions, and then the k-space will not be continuous, we obtain a discrete, densely spaced series of k vectors. The kdependent eigenenergies E n (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 F(t) 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 The second widely used option is the length (l) gauge, in which A l = 0, Φ l = −rF(t), obtained from the v gauge by using the gauge function Λ = −rA v (t). This yields (Note that both of these choices belong to the class of Coulomb gauges.) From now on, we shall use the notation A(t) ≡ A v (t). It should be kept in mind that electron wavefunctions in these gauges are connected by a space-time dependent unitary transfomation: (5) In the following we will consider pulsed excitation of duration T, and assume that both F(t) and 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 (r, t) is the identity]. Note that this choice avoids gauge-related ambiguities for t > T [18][19][20].
Under the effect of H g ext , any initial state will evolve in time, and can be expanded as a time dependent linear combination of the unperturbed eigenstates: As H g ext is gauge dependent, the expansion coefficients as well as the states |Ψ g (t) shall depend on the gauge chosen. We use the following notations in the two gauges of interest here: The following (invertible) relation holds between the expansion coefficients: 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 : 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 (k)] is practically unity for the valence band, and zero for the conduction bands. That is, to a very good approximation, we can write ρ(t = 0) = k |n 0 , k n 0 , 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 ρ l (t) = U (t)ρ v (t)U † (t) and the time evolution of the density operator is given by the von-Neumann equation . For the sake of simplicity, let us use the velocity gauge, when (in dipole approximation) the A 2 term in the Hamiltonian commutes with ρ v , and the remaining part does not mix states with different indices k (see the Appendix for more details). Therefore we have: where n |c nk (t)| 2 = 1, and the time dependence of these coefficients is determined by the von-Neumann equation [with the initial conditions of c nk (0) = δ n,n0 ]. Since the trace of the density operator is preserved during the dynamics, our choice of normalization means that Trρ equals the number 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. [21]), 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: 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 p kin , due to its direct connection to the velocity operator [22]. 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 In the length gauge A l =0, and we have therefore: 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 as the intraband contribution, while the microscopic, "polarization-like" parṫ as the interband contribution. Indeed, n, k|p kin |n, k /m is the band velocity in the nth band at point k, thus j (g) means the sum of all these velocities weighted by the corresponding populations. For comparison with textbook methods [23], let us note that n, k|p|n, k can also be calculated as m/ ∇ k E n (k), but Eqs. (11) and (12) treat inter-and intraband currents on an equal footing. Clearly, J = j g +Ṗ g , and it is gauge independent, but it does not necessarily hold for j g andṖ g separately. Let us investigate the relation (7) between the expansion coefficients c nk (t) and b nk (t). Assume that the dynamics is solved in the velocity gauge [i.e., c nk (t) are known functions of time], and let us transform the result to length gauge. We obtain: where k(t) = eA(t)/ − k ′ . The first consequence of this result is that whenever the velocity gauge density operator is diagonal in 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 k values being shifted by eA(t)/ . [This is in agreement with the acceleration theorem (1).] On the other hand, since the lattice-periodic functions u n,k are orthogonal only for the same index k, the integral in the last line of Eq. (15) is not proportional to δ nn ′ . This effect is closely related to the fact that the separation of J as the sum of intraband (j) and interband (Ṗ) components is not unambiguous. E.g., for j we have where N = Trρ. On the other hand, The last term in Eq. (16) is determined solely by the time dependence of the vector potential, while the coefficients c nk and b nk 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 j l (t) and 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 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 [3]. Having determined the Bloch-sates, the time dependence of j g andṖ g [the onedimensional 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 gaugeindependent 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" [24]. However, for optical frequencies and laser field strengths around 1 GV/m, the phase factor appearing in exp(−ieAr/ ) is almost constant over the size of an atom, its change is around π/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., A = 0), the separation J = j+Ṗ is clearly unambiguous, and e.g., 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. 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: where R denotes lattice vectors. In view of this, the matrix elements of H ext in velocity gauge [recall that A(t) ≡ A v (t)] are given by: 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 l ext (t) as given by −e 0 n, k|r|n ′ , k ′ F(t) cannot be simplified as above, because the operator r is not lattice-periodic.

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 2 2m −i ∂ ∂x + k 2 + V (x) u n,k (r) = E n (k)u n,k (r). (20) As numerical examples, we investigated two model potentials.
V 1 contains two localized attracting centers: , 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 [−π/2, π/2]. The second model potential is a modification of the one that has been used in Ref. [13]: , 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 {u n,k } N 1 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 = −π/a + (m − 1)∆k, m = 1, . . . , M and (M − 1)∆k = 2π/a. In this way, we obtain N × M functions u n,k , and -besides the matrix elements (19) -we can calculate their overlaps numerically. Note that -by construction -we have S kk nn ′ = δ nn ′ . 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 k(t) = eA(t)/ − 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 provided t ∈ [0, T ], and zero otherwise. Similarly to Ref. [3], we consider many-cycle, mid-infrared excitation. That is, λ = 3µ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.