A Gauge Invariant Formulation of Interband and Intraband Currents in Solids
AA Gauge Invariant Formulation of Interband and Intraband Currents in Solids
Guilmot Ernotte, T. J. Hammond, and Marco Taucer ∗ Joint Attosecond Science Laboratory, National Research Council ofCanada and University of Ottawa, Ottawa, Ontario K1A 0R6, Canada Department of Physics, University of Windsor, Windsor, Ontario N9B 3P4, Canada (Dated: August 29, 2018)Experiments and simulations in solid-state high harmonic generation often make use of the dis-tinction between interband and intraband currents. These two contributions to the total currenthave been associated with qualitatively different processes, as well as physically measurable signa-tures, for example in the spectral phase of harmonic emission. However, it was recently argued [P.F¨oldi, Phys. Rev. B , 035112 (2017)] that these quantities can depend on the gauge employedin calculations. Since physical quantities are expected to have gauge-independent values, this raisesthe question of whether the decomposition of the total current into interband and intraband con-tributions is physically meaningful, or merely a feature of a particular mathematical representationof nature. In this article, we explore this apparent ambiguity. We show that a closely related issuearises when calculating instantaneous band populations. In both the case of inter/intraband currentsand in the case of instantaneous band populations, we propose definitions which are gauge-invariant,and thus allow these quantities to be calculated consistently in any gauge. I. INTRODUCTION
The strong field of a pulsed laser can drive extremelynonlinear currents in a solid, leading to the emission ofhigh-order harmonics of the fundamental frequency thatcan span the visible spectrum and extend into the ex-treme ultraviolet [1]. The process can take place in a widerange of materials, from dielectrics to semiconductors tosemimetals, and can leave the material undamaged [2–7].Strong field and attosecond science in condensed matteris an extension of the long-standing field of high har-monic generation, which was for many years confined togas-phase atoms and molecules [8]. Solids remain a newarea of this field, in which some basic questions remainunanswered, while others may not even be precisely de-fined.In that respect, as part of the search for an underlyingphysical picture of the harmonic generation process, ex-perimentalists and theorists alike have focused much at-tention on the division of the total current into interbandand intraband processes. This conceptual separation isappealing in part because the interband picture bears astrong similarity to the well-understood gas-phase model[7, 9, 10], while the intraband picture is qualitatively dif-ferent and for the most part unique to the solid state[6, 11]. Experiments have access to the complex ampli-tude of the emitted harmonics which reflects the coherentsum of the interband and intraband contributions, what-ever their relative weight may be. For now, a clean sep-aration is only possible in calculations; any conclusionsabout the dominance of one mechanism or the other re-lies on a comparison with theoretical predictions of theinterband and intraband spectra [2, 6, 7, 10–15]. Re-cently, however, F¨oldi showed that this separation may ∗ [email protected] be gauge-dependent [16]. That is, a different choice ofthe gauge, which should leave all physical quantities un-changed, leads to different values for the interband andintraband currents. However, the total current, whichrelates to the experimentally observed harmonic spec-trum, is not gauge-dependent. This raises the questionof whether this conceptual decomposition of the currentis physically meaningful. Is the interband current an ob-servable?Other quantities may also be easy to calculate, buthard to access in experiments. An example is the in-stantaneous band population. In the strong-field physicsof solids, simulations often show a transient conductionband population which oscillates with the applied field,but mostly returns to the valence band at the end of thepulse. The fraction of the population which remains inthe conduction band at the end of the pulse depends onthe band structure and the pulse shape, as well as the de-phasing time constant [17]. While this final populationis gauge-independent, the transient population dynam-ics during the laser pulse’s illumination can be subjectto a gauge-dependence that is analogous to that of theinter- and intra-band currents, as we will show below.This raises the question of the significance of the instan-taneous band populations. Can such a quantity be pre-cisely defined, particularly given that in the presence ofa strong driving field the instantaneous eigenstates arethe dressed states, which are not the same as the field-free eigenstates? The question is all the more compellinggiven that some attosecond-probe experiments appear tomeasure precisely this quantity [18].The aim of the present article is to create gauge-invariant definitions of these quantities of interest. Westart by identifying the Hermitian operators correspond-ing to the instantaneous band populations, and to theinterband and intraband currents. Once defined, we de-rive their gauge-dependent transformations, which en-sures gauge-invariant physical predictions. While we pri- a r X i v : . [ c ond - m a t . m e s - h a ll ] A ug marily focus on the commonly used velocity and lengthgauges, our definition is equally valid in any other gauge.This paper is organized as follows. In Section II, weintroduce the theoretical formalism and the details of ournumerical calculations. Section III discusses an intuitive,but problematic, approach to defining interband and in-traband currents as well as band populations. In SectionIV, we provide a more rigorous definition of the instan-taneous band population, as a Hermitian operator. Wethen calculate the gauge-transformation of its matrix el-ements. Section V similarly describes the interband andintraband currents in terms of Hermitian operators, withcorresponding gauge transformations. The improved def-initions of these quantities yield gauge-invariant predic-tions. Finally, in Section VI, we show that our definitionsgive reasonable physical descriptions, and we discuss thechoices made in coming to this formulation. II. THEORETICAL APPROACH
For a single particle in a one-dimensional periodic po-tential, V ( x + a ) = V ( x ) with a lattice constant of a ,the Hamiltonian in the absence of the laser field isˆ H = ˆ p V (ˆ x ) . (1)Here and throughout this paper we use atomic units, ex-cept where other units are specified. The eigenstates ofthis Hamiltonian can be labeled by a band index, n , andthe crystal momentum, k :ˆ H | φ nk (cid:105) = ε n ( k ) | φ nk (cid:105) . (2)The energies, ε n ( k ), trace out the band structure, andthe Bloch functions, expressed in the position basis, havethe property (cid:104) x | φ nk (cid:105) ≡ φ nk ( x ) = (cid:112) a π e ikx u nk ( x ), where u nk ( x + a ) = u nk ( x ) is periodic and normalized over oneunit cell.As a model system, we use the previously studiedMathieu potential, V ( x ) = − V [1 + cos(2 πx/a )], with V = 0 .
37 and a = 8 atomic units [19–23]. We solvethe Time-Independent Schrodinger Equation (TISE) inthe position basis, with periodic boundary conditions, tofind the field-free eigenstates (Bloch states). These arethen used as the basis for calculations of the time dy-namics in a driving laser field. Figure 1 shows the bandstructure as a function of crystal momentum for the firstfive bands. The black circle in the center of the band withindex n = 1 represents the initial condition for simula-tions: | φ ,k =0 (cid:105) , a single electron at the Γ-point in band-1.Roughly speaking, bands 1 and 2 can be thought of as thevalence band and the first conduction band, respectively.However, our calculation considers band-0 to be unoccu-pied, as well as all other k -points in band-1. While thissimplification does not represent the reality of valencebands, it allows a comparison with previous reports andhas no effect on the conclusions of this work. -10020304010 1.0-0.5 0.0 0.5-1.0 n =0 n =1 n =2 n =3 n =4 FIG. 1. Band structure for the Mathieu potential, showingthe lowest five bands, labeled with their band indices, as afunction of crystal momentum in the first Brillouin zone. Theblack circle in band 1 shows the initial condition for our simu-lation. Orange and blue arrows/circles illustrate the dynamicsas described in the velocity and length gauges, respectively.Both gauges describe the time-dependent wavefunction as asuperposition of Bloch states at a particular crystal momen-tum, but in the velocity gauge the crystal momentum is fixed,while in the length gauge it oscillates with the vector poten-tial.
The effect of a laser field, which we treat here withinthe dipole approximation, leaves freedom with respectto the gauge chosen, since the field can be divided non-uniquely between a scalar potential, Φ, and a vector po-tential, A . In the length gauge, the field is incorporatedexclusively through the scalar potential,Φ ( l ) ( x, t ) = − F ( t ) · xA ( l ) ( t ) = 0 , (3)where F ( t ) represents the electric field. In the velocitygauge, the situation is reversed,Φ ( v ) ( x, t ) = 0 A ( v ) ( t ) = − (cid:90) t −∞ F ( t (cid:48) ) dt (cid:48) . (4)The Hamiltonian, including the interaction with the lightfield, under the dipole approximation, can be written inany gauge, labelled with superscript ( g ), asˆ H ( g ) = 12 (cid:104) ˆ p + A ( g ) ( t ) (cid:105) + V (ˆ x ) − Φ ( g ) (ˆ x, t ) . (5)Note that when the field and the (velocity gauge) vectorpotential are both zero (that is, before or after the pulse),the Hamiltonian reduces to ˆ H in both the velocity andlength gauge.Hermitian operators corresponding to observablestransform between the two gauges according to [24–26],ˆ O ( v ) = e − iA ( v ) ˆ x ˆ O ( l ) e iA ( v ) ˆ x , (6)and the wavefunctions are related by | ψ ( v ) (cid:105) = e − iA ( v ) ˆ x | ψ ( l ) (cid:105) . (7)More generally, wavefunctions and operators are trans-formed from gauge ( g ) to gauge ( g ) by the unitary op-erator ˆ U ( g ) → ( g ) ≡ e i [ A ( g − A ( g ]ˆ x . (8)In this article, we choose to equate the field-free opera-tor for our quantities of interest with the operator’s rep-resentation in the length gauge. While this is not theonly possible choice, we will provide a justification be-low, and show that this definition gives reasonable phys-ical predictions. Starting from the length gauge, then,the transformation to any other gauge, ( g ), is describedby ˆ U ≡ ˆ U ( l ) → ( g ) = e iA ( g ) ˆ x . (9)In each gauge, the vector potential dictates a time-dependent transformation. As long as wavefunctions andHermitian operators transform via the unitary opera-tor of Eq. 9, expectation values remain unchanged bythe gauge transformation. In other words, the differentgauges all represent the same physics.We solve the Time-Dependent Schr¨odinger Equation(TDSE) in the field-free basis numerically. We writewavefunctions in gauge g as | ψ ( g ) ( t ) (cid:105) = (cid:88) n (cid:90) BZ dk c ( g ) nk ( t ) | φ nk (cid:105) . (10)The amplitudes, c ( g ) nk ( t ), depend on the gauge. For nu-merical simulations, the laser field is defined by A ( v ) ( t ) = A cos ( ω t/ n c ) cos( ω t ), where A = 0 . n c = 11 is the number of cycles in thepulse, and a fundamental frequency of ω = 2 πc/λ , with c the speed of light and λ = 3 . µ m the wavelength.Our parameters are identical to those used by Wu et al. [19, 20].In the presence of the laser field, the Hamiltonian ex-pressed in the field-free basis acquires off-diagonal ele-ments that couple the Bloch states. This coupling is de-termined by the matrix elements of the momentum oper-ator in the velocity gauge, or by the matrix elements ofthe position operator in the length gauge. The resultingset of coupled differential equations for the coefficients, c ( g ) nk ( t ), takes a different form in each gauge. In the ve-locity gauge, i ∂∂t c ( v ) nk = (cid:20) ε n ( k ) + 12 A ( v )2 (cid:21) c ( v ) nk + A ( v ) (cid:88) n (cid:48) p nn (cid:48) ( k ) c ( v ) n (cid:48) k . (11) The momentum operator only couples states with thesame k -value, (cid:104) φ nk | ˆ p | φ n (cid:48) k (cid:48) (cid:105) ≡ p nn (cid:48) ( k ) δ ( k − k (cid:48) ), meaningthat the initial crystal momentum, in the velocity gaugedescription, remains constant even under the influence ofthe laser field. In the length gauge, i ∂∂t c ( l ) nk = ε n ( k ) c ( l ) nk + iF ∂∂k c ( l ) nk + F (cid:88) n (cid:48) ξ nn (cid:48) ( k ) c ( l ) n (cid:48) k . (12)The position matrix elements likewise contain a partwhich mixes states of the same k -value, denoted ξ nn (cid:48) ( k ),which, for non-degenerate states, is given by ξ nn (cid:48) ( k ) = − ip nn (cid:48) ( k ) / ( ε n ( k ) − ε n (cid:48) ( k )). However, the position oper-ator additionally contains a differential term which cou-ples neighbouring k -values, and leads to the accelerationtheorem: a state initially having k = k evolves intostates with k ( t ) = k + A ( v ) ( t ). The acceleration theo-rem illustrates an important difference between these twogauges: in the velocity gauge, the crystal momentum re-mains fixed, while in the length gauge, the evolution ofthe wavefunction in the laser field leads to an oscillationof the electron’s crystal momentum. This is illustratedby the coloured circles and arrows in Fig. 1. III. DEFINITIONS BASED ON BAND INDICES
In this section we discuss an intuitive yet problem-atic procedure for defining band populations and in-ter/intraband currents, in which the band indices of co-efficients are used to identify band-dependent quantities.This procedure is known to give gauge-dependent results[16]. In subsequent sections, we will attempt to improveupon these definitions.Solving the Schr¨odinger equation in either gauge givesthe time-dependent coefficients of the basis states, andthereby the wavefunction. In particular, the coefficients, c ( g ) nk ( t ), are associated with the state | φ nk (cid:105) , and, accordingto the Born rule, their modulus squared seems to repre-sent the probability of finding an electron in that partic-ular state. This reasoning implies that the instantaneousband population in band n is (cid:82) dk | c ( g ) nk | . The result-ing instantaneous conduction band population is shownin Fig. 2a, for the velocity gauge (orange) and lengthgauge (blue). Both gauges describe a transient popula-tion which oscillates and nearly completely returns to thevalence band at the end of the pulse. Both calculationsagree on the final population that remains in the con-duction band at the end of the pulse. However, in thevelocity gauge, the transient conduction band populationis much larger, as much as 40%, and it is peaked at the(velocity gauge) vector potential maxima, whereas in thelength gauge the apparent conduction band population ispeaked at the vector potential zeros. Also, whenever thevector potential (Fig. 2c) is zero, indicated by verticaldashed lines, the two gauges agree. This is to be expectedsince the unitary transformation (Eq. 9) is the identityin that case. But in general, the two calculations providevery different pictures of the conduction band popula-tion, and whenever the vector potential is non-zero, it isnot clear which one to trust. In the recent literature, itappears that both quantities have been reported [18, 27–30]. It would be desirable to find a gauge-independentformulation of the instantaneous band population.Turning now to the current, we can express it simplyin terms of the kinematical momentum, and expand itsexpectation value in terms of Bloch states, j ( t ) = −(cid:104) ψ ( g ) ( t ) | ˆ p ( g ) kin | ψ ( g ) ( t ) (cid:105) = − (cid:88) n (cid:90) BZ dk | c ( g ) nk | (cid:104) φ nk | ˆ p ( g ) kin | φ nk (cid:105)− (cid:88) n,n (cid:48) (cid:54) = n (cid:90) BZ dkc ∗ ( g ) nk c ( g ) n (cid:48) k (cid:104) φ nk | ˆ p ( g ) kin | φ n (cid:48) k (cid:105) , (13)noting that in the velocity gauge the kinematical momen-tum is ˆ p ( g ) kin = ˆ p + A ( g ) , where ˆ p refers to the canonical mo-mentum. The summation over all basis states has beensplit into a summation over terms involving the sameband indices (intraband) and a summation over termsinvolving different band indices (interband). Figure 2bshows the interband current, so defined, as calculatedin the velocity (orange) and length (blue) gauges. Asin the case of the instantaneous conduction band pop-ulation, the two calculations give very different results.Importantly, the total current, shown in black, is gauge-independent as long as the kinematical momentum isused. Since the total current is what gives rise to the mea-sured harmonic spectrum, there is no question that it is aphysically meaningful quantity, and its gauge-invarianceis expected. However, as long as the treatment of inter-band and intraband currents depends on the gauge, it isnot clear that these are physical quantities. IV. BAND POPULATIONS
In this section we propose a gauge-invariant definitionof the instantaneous band populations. Our approachis to identify its corresponding Hermitian operator. Wethen transform it with the proper unitary operators tomake it gauge-independent. Before proceeding, we con-sider the general question of the population of any quan-tum state, | S (cid:105) . The probability to find a system in thisstate is evidently the squared modulus of the projectionof the system’s state onto | S (cid:105) . Another way to say thisis that it is the expectation value of a Hermitian op-erator (an observable), namely the projection operator,ˆΠ S ≡ | S (cid:105)(cid:104) S | .Likewise the operator representing band population isan operator that projects onto all the states within agiven band, m ,ˆΠ m = (cid:90) BZ dq ˆΠ mq ≡ (cid:90) BZ dq | φ mq (cid:105)(cid:104) φ mq | , (14) P opu l a t i on I n t e r band C u rr en t A ( t ) Time (fs)VelocityLength x10 abc
FIG. 2. Dynamics calculated using definitions based on bandindices. (a) Conduction band population, defined as thesquared modulus of the time-dependent coefficients, calcu-lated in the velocity and length gauges. The region on theright side is magnified in the vertical dimension by a factor often. (b) Interband current, defined as the second summationin Eqn. 13, calculated in the velocity gauge and length gauge.The total current, which is gauge-independent provided onemakes the transformation ˆ p ( g ) = ˆ p + A ( g ) , is shown in blackfor reference. (c) Velocity gauge vector potential as a functionof time. Vertical dashed lines in all plots indicate the zeros ofthis vector potential. where Π mq is the projection onto a single eigenstate, | φ mq (cid:105) , of the field-free Hamiltonian. Since this opera-tor was defined without reference to a gauge, nor a field,we refer to it as a “field-free operator”. The expectationvalue of ˆΠ m represents the band population, and Eq. 9provides its transformation. In a gauge ( g ) the operatortransforms as ˆΠ ( g ) mq = e − iA ( g ) ˆ x ˆΠ mq e iA ( g ) ˆ x , (15)and its matrix elements, needed to compute the expec-tation value, in the field-free basis are (cid:104) φ nk | ˆΠ ( g ) mq | φ n (cid:48) k (cid:48) (cid:105) = (cid:104) φ nk | e − iA ( g ) ˆ x | φ mq (cid:105)(cid:104) φ mq | e iA ( g ) ˆ x | φ n (cid:48) k (cid:48) (cid:105) = U † mn ( q, k ) U mn (cid:48) ( q, k (cid:48) ) . (16)The matrix elements of the transformation operatorcan be shown to be U mn ( q, k ) = δ ( k + A ( g ) − q )∆ mn ( k + A ( g ) , k ) , (17)where we have introduced the function∆ nm ( k , k ) ≡ (cid:104) u nk | u mk (cid:105) , (18) Band 0Band 1Band 2 b B and I nde x a FIG. 3. Squared modulus of the functions, ∆ n, (0 , k ), plottedover ten Brillouin zones of reciprocal space, for (a) the lowestseven bands and (b) for the lowest three bands. These indi-cate the degree to which each periodic function, at k = 0, isrequired to represent a wavefunction with crystal momentum k in band 1. which is the overlap integral of the periodic parts of theBloch wavefunctions within one unit cell. While the totalwavefunctions, | φ nk (cid:105) , are all mutually orthogonal, this isnot true for the periodic parts, | u nk (cid:105) . Whenever k = k ,the function ∆ nm ( k , k ) reduces to a Kronecker deltafunction, δ nm , but not otherwise. This is a statementof the fact that, at a particular value of the crystal mo-mentum, k = k , the periodic functions form a com-plete orthonormal set, u nk ( x ), within the Hilbert spaceof a single unit cell. But a different choice of k leads toa different set of functions, which will be mutually or-thogonal, but need not be orthogonal to the functions ofthe first set. This can be seen in Fig. 3, which shows | ∆ n, (0 , k ) | , as a function of k for different bands, n .At k = 0, ∆ (0 ,
0) = 1 while all other matrix elementsare zero. However, as k departs from zero, neighbouringbands become important (Fig. 3b), and as k increases(over several Brillouin zones), higher-lying bands havethe largest contributions (Fig. 3a). At a glance, thisgives a sense of how many bands are required at the Γ-point to express the periodic part of the wavefunction ata different point, k , in reciprocal space. We note thatthese ∆ functions have been discussed previously in thesolid-state literature [31].From the Dirac delta function of Eq. 17, we imme-diately see that in the velocity gauge the operator that reports the population of the state with crystal momen-tum q in band m involves matrix elements with a differentcrystal momentum, q − A ( g ) , and may involve matrix ele-ments in all bands. This is why band populations appear to be very different when calculated in the two gauges.Finally, the matrix elements of the band projectionoperator, found by integrating the wavevector q over oneBrillouin zone, are (cid:104) ˆΠ ( g ) m (cid:105) nn (cid:48) ,kk (cid:48) = ∆ ∗ mn ( k + A ( g ) , k )∆ mn (cid:48) ( k + A ( g ) , k ) δ ( k − k (cid:48) ) . (19)The expectation value of this projection operator rep-resents our proposed formulation of instantaneous bandpopulations. Figure 4a shows the conduction band pop-ulation calculated using Eq. 19 in the velocity gauge(dashed orange) and length gauge (filled blue). The twocalculations overlap exactly, showing that this definitionindeed gives a gauge-invariant value of the instantaneousband population.The conduction band population exhibits sharp peaksat the vector potential zeros (near the peaks of the elec-tric field), shown more clearly in Fig. 4b. The sharpnessof these peaks is in part due to the fact that the valenceband is only occupied at a single value of k . A filledvalence band will lead to broader transient peaks in theconduction band population. The maximum transientpopulation transfer is about 5%, much lower than the40% described by the previous velocity-gauge calculationof Fig. 2a. We also note that, since the transformationoperator for the length gauge is the identity operator, thelength gauge calculations of Figs. 4a and 2a are identical.We will revisit this fact in Section VI. V. INTERBAND AND INTRABANDCURRENTS
Having considered the band populations during thepulse, and the closely related projection operators, weare in a position to revisit the question of interbandand intraband currents. As before, we start by identify-ing a Hermitian operator, and then determine its gauge-dependent transformation. To do this, we again make useof projection operators to define the field-free operatorsˆ j ra = − (cid:88) n ˆΠ n ˆ p kin ˆΠ n (20)and ˆ j er = − (cid:88) n,n (cid:48) (cid:54) = n ˆΠ n ˆ p kin ˆΠ n (cid:48) . (21)The expectation values of these operators, in the absenceof any applied fields, reproduce precisely the two summa-tions shown in Eq. 13. However, when a field is applied, agauge must be chosen and the operators transformed ap-propriately. The transformation to the length gauge is,again, trivial and leaves the decomposition unchanged. P opu l a t i on a VelocityLength d -5 0 5-0.10.00.1 I n t e r band C u rr en t Time (fs) c b -5 0 50.040.020.00 FIG. 4. (a) Conduction band population, defined using Eq.19, as a function of time for velocity and length gauge, whichoverlap. (b) Magnified view of the conduction band popula-tion, showing transient peaks at the vector potential zeros,as well as high-frequency oscillations. (c) Interband current,defined using Eq. 25, as a function of time for velocity andlength gauge. (d) Magnified view of interband current showsthat the two gauges agee within numerical acuracy.
However, when transforming the operator to the veloc-ity gauge, we must transform not only the momentumoperator, but also the projection operators. By simplydividing the current according to the terms shown in Eq.13, we perform the required transformations on the mo-mentum operator (ˆ p ( g ) = ˆ p + A ( g ) ) and the wavefunctions(Eq. 7), but we fail to take into account the transforma-tion of the projection operators.The gauge-transformed intraband current operator isˆ j ( g ) ra = − (cid:88) m (cid:90) (cid:90) BZ dqdq (cid:48) ˆΠ ( g ) mq (ˆ p + A ( g ) ) ˆΠ ( g ) mq (cid:48) . (22)Its matrix elements are (cid:104) ˆ j ( g ) ra (cid:105) nn (cid:48) ,kk (cid:48) = − (cid:88) mll (cid:48) ∆ ∗ mn ∆ ml ∆ ∗ ml (cid:48) ∆ mn (cid:48) × (cid:104) p ll (cid:48) ( k ) + δ ll (cid:48) A ( g ) (cid:105) δ ( k − k (cid:48) ) , (23)where the arguments of the functions, ∆, are understoodto be ∆ ab ≡ ∆ ab ( k + A ( g ) , k ).Likewise, for the interband current, we haveˆ j ( g ) er = − (cid:88) m,m (cid:48) (cid:54) = m (cid:90) (cid:90) BZ dqdq (cid:48) ˆΠ ( g ) mq (ˆ p + A ( g ) ) ˆΠ ( g ) m (cid:48) q (cid:48) , (24) whose matrix elements are (cid:104) ˆ j ( g ) er (cid:105) nn (cid:48) ,kk (cid:48) = − (cid:88) mm (cid:48) ll (cid:48) ∆ ∗ mn ∆ ml ∆ ∗ m (cid:48) l (cid:48) ∆ m (cid:48) n (cid:48) × (cid:104) p ll (cid:48) ( k ) + δ ll (cid:48) A ( g ) (cid:105) δ ( k − k (cid:48) ) , (25)where the summation omits terms with m = m (cid:48) . Fig-ures 4c and d show the interband current calculated usingEquation 25 in the velocity gauge (dashed orange) andlength gauge (solid blue). Here again, the agreement be-tween the two calculations shows that our definition isgauge-invariant. VI. DISCUSSION
We have proposed new formulations for the interbandand intraband currents, as well as the instantaneous bandpopulations. This was done by identifying Hermitian op-erators corresponding to each of these, and then apply-ing a gauge-dependent unitary transformation. We haveshown that the resulting operators give gauge-invariantpredictions for these quantities, providing a possible res-olution to the issue of gauge-dependence. In addition tothis, our formulation allows all of these quantities to becomputed from any gauge. In particular, the velocitygauge may present advantages regarding computation insome cases [24].The harmonic spectrum for the strongly driven Math-ieu potential, calculated entirely in the velocity gauge,is shown in Fig. 4. The total current (black line) dis-plays the previously discussed double plateau. The dif-ficulty in retrieving interband and intraband informa-tion was noted by Wu et al. [19], who suggested anapproach involving projection onto Houston states. Weagree with their approach and conclusions, however, sincethey project onto a time-dependent basis, they effectivelychange gauges and derive a new equation of motion forthe wavefunction. We discuss this further in AppendixA. Our approach allows us to consider the separationinto interband (not shown) and intraband (red line) con-tributions, entirely in the velocity gauge. Furthermore,it is possible to make finer divisions of the current. Forinstance, we may define the interband current betweenbands one and two as ˆ j = ˆΠ ˆ p kin ˆΠ + h . c . , whose spec-trum is shown as the blue shaded region. This currentaccounts for the first plateau. Likewise, the summed in-terband currents from band one to bands three and four,ˆ j + ˆ j , shown as the pink shaded region, account forthe second plateau. The relation of these plateaus tointerband dynamics was previously argued based on theenergy ranges at play and the intensity scaling [20, 21],but here we can compute the interband dynamics directlyeven though our calculation is in the velocity gauge.It is worth noting that the gauge-invariant quantitiesthat we defined coincide with the earlier problematic def-initions in the length gauge (those of Section III). In TotalIntraband j j + j -1 -4 -7 -10 Harmonic Order S pe c t r u m FIG. 5. Harmonic spectrum for the strongly driven Mathieupotential. The black line shows the spectrum of the totalcurrent. The red line shows that of the intraband currentwhich makes a negligible contribution for harmonics above theseventh order. The blue shaded region shows the spectrum ofthe interband current due to transitions between bands oneand two. Likewise, the pink shaded region shows the spectrumof the combined interband currents due to transitions betweenbands one and thre, and bands one and four. All quantitiesshown are calculated in the velocity gauge. our definitions, the field-free operators are transformedby the unitary operator e iA ( g ) x . However, in the lengthgauge, since A ( l ) = 0, our approach leaves the field-freeoperator unchanged. But could a different choice havebeen made? In particular, could we not have chosena definition such that the velocity gauge operator cor-responds to the field-free operator? Yes, such a choicecould have been made, and the transformations of Eq. 8would ensure gauge invariance. The resulting populationoperator in gauge ( g ) would be e − i [ A ( g ) − A ( v ) ] x ˆΠ mq e i [ A ( g ) − A ( v ) ] x , (26)to be compared with our proposed definition, Eq. 15.What, then, justifies the choice to associate the lengthgauge operator with the field-free operator? While weare not aware of any hard and fast justification for ourchoice, we suggest two possible arguments here, the firstesthetic and the second an appeal to reasonableness.Esthetically, the alternate definition of Eq. 26 de-fines the operator in gauge ( g ) by making reference tothe vector potential in that gauge and in gauge ( v ), aswell as referencing the field free operator. Thus, in anygauge, this alternate definition must always refer backto the velocity gauge. In contrast, our proposed defini-tion, Eq. 15, describes the operator in gauge ( g ) only interms of the field free operator and the vector potentialin that gauge. Perhaps more importantly, our definitions give reasonable answers, matching the understanding de-scribed in the literature. Figure 5 shows that our gauge-invariant formulation explains the two plateaus in theharmonic spectrum as arising from different interbandcurrents, consistent with previous arguments. A morestraightforward example is the response of the system toa constant electric field. For a sufficiently weak field, oneexpects an adiabatic evolution of the wavefunction withina single band, together with Bloch oscillations in the in-traband current. Appendix B shows that this behaviouris captured by our definitions, but not by alternative for-mulations like Eq. 26.The questions we addressed in this article attempt todefine instantaneous band-dependent quantities in thepresence of a driving field. Yet the bands themselvesare, in some sense, a field-free concept. A considera-tion of laser-dressed states is therefore an important partof this discussion. Laser-dressed states can be definedin various ways, either as instantaneous eigenstates ofa Hamiltonian, which will depend on the chosen gauge,in a cycle-averaged (Floquet) way, or otherwise [17, 32–34]. Here, we have not touched on these issues, asidefrom the Houston states, which are the instantaneouseigenstates of the velocity gauge Hamiltonian. In atomicand molecular systems, the dressed states have played animportant role in resolving unphysical anomalies in cal-culations, which are often gauge-dependent [33, 34]. Wemay then expect such considerations to shed further lighton questions of, for example, instantaneous state popu-lations in a laser field. We leave a detailed considerationof dressed states, as they relate to gauge-dependence, forfuture study.To conclude, we have proposed a formulation of instan-taneous band populations, and interband and intrabandcurrents which gives gauge invariant predictions. Our ap-proach was to define Hermitian operators correspondingto these quantities, which can be done without makingreference to a field or a gauge. These field-free operatorsare then transformed by the well-known unitary opera-tor that connects gauges and ensures gauge invariance.We have demonstrated numerically that this gives iden-tical results in the velocity and length gauges, and thatthese results are consistent with the community’s use ofthese terms. The results give these important quanti-ties a more rigorous definition, which, by removing thegauge-dependence, establishes that these are physicallymeaningful. ACKNOWLEDGEMENTS
For helpful discussions and comments, we would like tothank our colleagues, Michael Spanner, Thomas Brabec,Chris McDonald, Andr´e Staudte, and Paul B. Corkum.This work was supported by the Vanier Canada GraduateScholarship program (G.E.) and by the Air Force Office ofScientific Research Multidisciplinary University ResearchInitiative grant number FA9550-15-1-0037.
Appendix A: Relation to Houston States
Houston states have played an important role in thedevelopment of solid-state physics. They were originallyproposed as approximate solutions of the length gaugeHamiltonian in a constant electric field [35]. Later, theywere employed by Krieger and Iafrate to find an analyt-ical formulation of the TDSE in the velocity gauge [36].This latter approach, starting from the velocity gauge,was used by Wu et al. to make the separation betweeninterband and intraband currents [19]. They showed thatwhen the TDSE in the velocity gauge is solved in the ba-sis of Houston states, the simple decomposition accordingto band indices gives sensible results. Indeed, it has longbeen noted that there is a close connection between theHouston state basis, and gauge transformations. In thisAppendix, we address this connection.In velocity gauge calculations, the Houston states aredefined as | (cid:101) φ nk (cid:105) = e − iA ( v ) x | φ nk ( t ) (cid:105) , (A1)where k ( t ) = k + A ( v ) , and the vector potential is un-derstood to be that of the velocity gauge. We make twoimportant remarks on these states. First, we note thatthe crystal momentum (on the right side of the equation)acquires a time dependence which is a consequence ofthe acceleration theorem. Thus, the states | φ nk ( t ) (cid:105) canbe thought of as “accelerated Bloch states”. They areeigenstates of H , albeit with a time-dependent crystalmomentum. Second, the the phase factor which mul-tiplies this state is precisely the unitary transformationthat describes a gauge transformation, U † . Thus, theHouston states defined by Eq. A1 can be thought ofas the velocity gauge representation of accelerated Blochstates. Indeed, the original treatment by Houston didnot include the transformation operator, e − iA ( v ) x , sinceit employed the length gauge.Considering the TDSE in the Houston-state basis, nu-merically, one solves for the complex amplitudes of thebasis states: (cid:104) (cid:101) φ nk | ψ ( v ) (cid:105) = (cid:110) (cid:104) φ nk ( t ) | e iA ( v ) x (cid:111) · (cid:110) e − iA ( v ) x | ψ ( l ) (cid:105) (cid:111) = (cid:104) φ nk ( t ) | ψ ( l ) (cid:105) . (A2)That is to say the coefficients describing the Houstonstates in the velocity gauge are identical to those de-scribing accelerated Bloch states in the length gauge. Itis thus not surprising that the coupled differential equa-tions describing these coefficients end up being essentiallyidentical to the length gauge Schr¨odinger equation. In-deed, one arrives at the same differential equations in thelength gauge by using an accelerated frame [12].The approach of using a Houston-state basis does allowa decomposition into interband and intraband currents,however it involves a change in the differential equationsemployed, effectively requiring a return to the lengthgauge. In this sense, it does not directly address the issue of gauge-dependence. Also, it does not take advantageof the potential computational benefits of the velocitygauge. Our approach allows the TDSE to be solved en-tirely in the velocity gauge using the field-free basis, andstill allows an unambiguous determination of interbandand intraband currents. Furthermore, by defining quan-tities as Hermitian operators, we provide a more rigorousformulation which is equally valid in any gauge. Appendix B: Description of Bloch Oscillations inVelocity Gauge
The gauge-invariant definitions that we have proposedare chosen by first identifying the Hermitian operatorscorresponding to the observables of interest. This doesnot require us to think about a gauge, or even a field forthat matter: currents or band populations are quantitiesthat can be associated with the state of a system evenin the absence of a field. Once these Hermitian opera-tors are defined, gauge-invariance comes from the unitarytransformation of Equation 9. This is the thought pro-cess that leads to the definitions in Eqs. 15, 22, and 24.However, as noted above, the choice of a gauge-invariantdefinition based on a field-free operator is not unique.In this appendix, we use the simple case of Bloch oscil-lations to argue that our proposed definitions coincidewith the commonly accepted and discussed behavioursof band populations and currents.The definitions that we propose cause the length-gaugeoperator to be equal to the field free operator. An alter-nate definition for the instantaneous band population isdescribed by Eq. 26, which instead causes the velocitygauge operator to coincide with the field-free operator.One can likewise put forward an alternate definition forthe intraband current, e − i ( A ( g ) − A ( v ) )ˆ x ˆ j ra e i ( A ( g ) − A ( v ) )ˆ x , (B1)and for the interband current, e − i ( A ( g ) − A ( v ) )ˆ x ˆ j er e i ( A ( g ) − A ( v ) )ˆ x . (B2)All of these alternate definitions also give gauge-invariantanswers, like the ones we propose. However, the dynam-ics they describe do not match common expectations foran electron in a constant electric field. In a constant elec-tric field, one expects an electron to evolve adiabaticallywithin a single band (for sufficiently low field), undergo-ing Bloch oscillations.Figure 6 shows the dynamics of the Mathieu potentialin a constant field, comparing our proposed definitions tothe alternate definitions of Eqs. 26, B1, and B2. The con-stant field of − . × − turns on exponentially, witha time constant of 4 . k = 0. Time (fs) P opu l a t i on C u rr en t
80 120400.00.51.0-4-2024 0
Band 1ProposedAlternate Band 2 Band 3intra interAlternateProposed ab FIG. 6. Dynamics of the Mathieu potential in a static elec-tric field, comparing our proposed definitions (in blue) withthe alternate definitions of Eqs. 26, B1, and B2 (in orange).(a) Interband and intraband currents. Bloch oscillations canbe seen, and the Bloch period is shown. (b) Band populationsas a function of time. In the proposed definition, the popu-lation remains in band 1 throughout, reflecting the adiabaticevolution of the state.
In Fig. 6a, the proposed definition of the intraband(dashed blue) clearly exhibits Bloch oscillations, whilethe interband (solid blue) current is zero, as expected.The current described by the alternate definitions (in or-ange) is noticeably different. Most importantly, the al-ternate definitions describe Bloch oscillations as an inter-band current, contrary to the community’s understand-ing. The band populations described by our proposeddefinitions also match the common understanding. InFig. 6b, we see that the proposed definitions (blue) de-scribe a population which remains entirely in band 1,with the other band populations remaining zero. Thisdescribes the adiabatic evolution of the electron within asingle band. On the other hand, the alternate definitionshows that the population moves from band 1 to bands 2and 3, and eventually to higher-lying bands (not shown).The numerical results of Fig. 6 show that our proposeddefinitions give the correct results in a simple case wherethere is a consensus on the expected physical behaviour. [1] S. Ghimire, A. D. DiChiara, E. Sistrunk, P. Agostini,L. F. DiMauro, and D. A. Reis, Nature Physics , 138(2011).[2] O. Schubert, M. Hohenleutner, F. Langer, B. Urbanek,C. Lange, U. Huttner, D. Golde, T. Meier, M. Kira, S. W.Koch, and R. Huber, Nature Photonics , 119 (2014).[3] H. Liu, Y. Li, Y. S. You, S. Ghimire, T. F. Heinz, andD. A. Reis, Nature Physics , 262 (2017).[4] M. Sivis, M. Taucer, G. Vampa, K. Johnston, A. Staudte,A. Yu. Naumov, D. M. Villeneuve, C. Ropers, and P. B.Corkum, Science , 303 (2017).[5] N. Yoshikawa, T. Tamaya, and K. Tanaka, Science ,736 (2017).[6] T. T. Luu, M. Garg, S. Yu. Kruchinin, A. Moulet, M. Th.Hassan, and E. Goulielmakis, Nature , 498 (2015).[7] G. Vampa, T. J. Hammond, N. Thir´e, B. E. Schmidt,F. L´egar´e, C. R. McDonald, T. Brabec, and P. B.Corkum, Nature , 462 (2015).[8] F. Krausz and M. Ivanov, Reviews of Modern Physics , 163 (2009).[9] P. B. Corkum, Physical Review Letters , 1994 (1993).[10] G. Vampa, C. R. McDonald, G. Orlando, P. B. Corkum,and T. Brabec, Physical Review B , 064302 (2015).[11] M. Garg, M. Zhan, T. T. Luu, H. Lakhotia, T. Kloster-mann, A. Guggenmos, and E. Goulielmakis, Nature ,359 (2016).[12] G. Vampa, C. R. McDonald, G. Orlando, D. D. Klug, P. B. Corkum, and T. Brabec, Physical Review Letters , 073901 (2014).[13] T. J. Hammond, S. Monchoc´e, C. Zhang, G. Vampa,D. Klug, A. Yu. Naumov, D. M. Villeneuve, and P. B.Corkum, Nature Photonics , 594 (2017).[14] M. Hohenleutner, F. Langer, O. Schubert, M. Knorr,U. Huttner, S. W. Koch, M. Kira, and R. Huber, Nature , 572 (2015).[15] M. Taucer, T. J. Hammond, P. B. Corkum, G. Vampa,C. Couture, N. Thir´e, B. E. Schmidt, F. L´egar´e, H. Selvi,N. Unsuree, B. Hamilton, T. J. Echtermeyer, and M. A.Denecke, Physical Review B , 195420 (2017).[16] P. F¨oldi, Physical Review B , 035112 (2017).[17] C. R. McDonald, G. Vampa, P. B. Corkum, andT. Brabec, Physical Review Letters , 173601 (2017).[18] A. Sommer, E. M. Bothschafter, S. A. Sato, C. Jakubeit,T. Latka, O. Razskazovskaya, H. Fattahi, M. Jobst,M. Schweinberger, V. Shirvanyan, V. S. Yakovlev,R. Kienberger, K. Yabana, N. Karpowicz, M. Schultze,and F. Krausz, Nature , 86 (2016).[19] M. Wu, S. Ghimire, D. A. Reis, K. J. Schafer, and M. B.Gaarde, Physical Review A , 043839 (2015).[20] M. Wu, D. A. Browne, K. J. Schafer, and M. B. Gaarde,Physical Review A , 063403 (2016).[21] L. Liu, J. Zhao, W. Dong, J. Liu, Y. Huang, and Z. Zhao,Physical Review A , 053403 (2017).[22] T. Ikemachi, Y. Shinohara, T. Sato, J. Yumoto, M. Kuwata-Gonokami, and K. L. Ishikawa, Physical Re-view A , 043416 (2017).[23] X. Liu, X. Zhu, P. Lan, X. Zhang, D. Wang, Q. Zhang,and P. Lu, Physical Review A , 063419 (2017).[24] Y. C. Han and L. B. Madsen, Physical Review A ,063430 (2010).[25] A. D. Bandrauk, F. Fillion-Gourdeau, and E. Lorin,Journal of Physics B: Atomic, Molecular and OpticalPhysics , 153001 (2013).[26] We note that the Hamiltonian operator does not trans-form like the operators of other Hermitian observable. Itstransformation is found by substituting the transformedwavefunction into the TDSE.[27] K. S. Virk and J. E. Sipe, Physical Review B , 035213(2007).[28] N. Tancogne-Dejean, O. D. M¨ucke, F. X. K¨artner, andA. Rubio, Nature Communications , 745 (2017).[29] M. Schultze, E. M. Bothschafter, A. Sommer, S. Holzner, W. Schweinberger, M. Fiess, M. Hofstetter, R. Kien-berger, V. Apalkov, V. S. Yakovlev, M. I. Stockman, andF. Krausz, Nature , 75 (2013).[30] F. Schlaepfer, M. Lucchini, S. A. Sato, M. Volkov,L. Kasmi, N. Hartmann, A. Rubio, L. Gallmann, andU. Keller, Nature Physics , 560 (2018).[31] N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, andD. Vanderbilt, Reviews of Modern Physics , 1419(2012).[32] T. Higuchi, M. I. Stockman, and P. Hommelhoff, Phys-ical Review Letters , 213901 (2014).[33] O. Smirnova, M. Spanner, and M. Ivanov, Journal ofPhysics B: Atomic, Molecular and Optical Physics ,S307 (2006).[34] O. Smirnova, M. Spanner, and M. Ivanov, Journal ofModern Optics , 1019 (2007).[35] W. V. Houston, Physical Review , 184 (1940).[36] J. B. Krieger and G. J. Iafrate, Physical Review B33