Floquet-Theoretical Formulation and Analysis of High-Harmonic Generation in Solids
FFloquet-Theoretical Formulation and Analysis of High-Harmonic Generation in Solids
Tatsuhiko N. Ikeda, Koki Chinzei, and Hirokazu Tsunetsugu
Institute for Solid State Physics, University of Tokyo, Kashiwa, Chiba 277-8581, Japan (Dated: November 13, 2018)By using the Floquet eigenstates, we derive a formula to calculate the high-harmonic componentsof the electric current (HHC) in the setup where a monochromatic laser field is turned on at sometime. On the basis of this formulation, we study the HHC spectrum of electrons on a one-dimensionalchain with the staggered potential to study the effect of multiple sites in the unit cell such as thesystems with charge density wave (CDW) order. With the help of the solution for the Floqueteigenstates, we analytically show that two plateaus of different origins emerge in the HHC spectrum.The widths of these plateaus are both proportional to the field amplitude, but inversely proportionalto the laser frequency and its square, respectively. We also show numerically that multi-step plateausappear when both the field amplitude and the staggered potential are strong.
I. INTRODUCTION
High-harmonic generation (HHG) is the basis for theattosecond physics and has attracted renewed attentionowing to its successful observations in bulk solids drivenby a strong laser field [1–9]. These observations have re-vealed that the HHG in solids has characteristics differentfrom those in atomic gases [10]. For example, the high-energy cutoff of the output spectrum scales linearly withthe input field amplitude [1, 2, 4] rather than its squarei.e. the laser intensity. Besides, multiple plateaus emergein the high-harmonic output spectrum for very large am-plitude [5]. To understand the microscopic mechanism ofthese unique features of HHG in solids, many theoreticaland experimental studies are actively being conducted.A theoretical approach to this problem is to analyzethe electron dynamics in solids in the time domain. Inthis approach, the two-band semiconductor Bloch equa-tion [11–16] and the time-dependent Schr¨odinger equa-tion [17–23] were numerically solved in the presence ofa pulse electric field, and analysis was performed for thehigh-harmonic components of the electric current (HHC),which work as the source of the HHG. For both equa-tions, the unique scaling in solids is reproduced in one-dimensional models and it has been shown that the in-terband transition plays an important role as well as theintraband dynamics. Thus the time-domain approachhas successfully reproduced the experimental observa-tions, but it does not fit analytical approaches and it isnot straightforward to obtain systematic understandingof microscopic physics. For instance, it is obscure whythe output spectrum has peaks at multiples of the inputlaser frequency since the pulse input has a continuousspectrum.A complementary theoretical approach is to invoke theFloquet theory [24] and analyze the electron dynam-ics in the frequency domain. In this approach, the in-put electric field is idealized to have an exact period-icity in time, and this periodicity is utilized to definethe Floquet eigenstates [25], which correspond to the so-lutions of the time evolution equation. For the time-dependent Schr¨odinger equation, early studies [26–29]analyzed the Floquet eigenstates and the HHC carried by them. More recently, the Floquet theory has been ap-plied to one-dimensional systems [30–33], graphene andcarbon nanotubes [34–36], and three-dimensional sys-tems [37]. However, it has not been discussed well howthe Floquet eigenstates are related to the initial statesin recent experiments. In addition, those previous stud-ies are mostly numerical, and the characteristics of theHHG in solids have not been fully understood. There arealso Floquet-theoretical approaches for the semiconduc-tor equation [38, 39]. Higuchi et al. [38] considered thequasistatic limit of the input electric field and discussedthe HHC originating from the Bloch oscillation. The highharmonics induced by this mechanism are multiples ofthe Bloch frequency Ω B , not the input laser frequency Ω,and this regime differs from that of the experiments inRef. [1, 5].In this paper, we develop a theoretical frameworkbased on the Floquet theory to investigate the mecha-nism of the HHG in solids. Considering the setup thatan ac electric field with frequency Ω is turned on atsome time to drive the system, we derive a formula toobtain the HHC spectrum from the Floquet eigenstates[see Eqs. (25) and (26)]. On the basis of our formula-tion, we then analyze the HHC spectrum of electrons ona one-dimensional chain [see Eq. (1)] to study the effect ofmultiple sites in the unit cell. We realize a two-site unitcell by introducing a staggered potential, which changesits sign alternately along the chain. In the absence of thestaggered potential, we see the presence of a plateau inthe HHC spectrum for strong field. Then we analyticallyshow that the staggered potential induces another widerplateau, which sets in already for a weaker field. We ar-gue that the widths of both plateaus scale linearly withthe field amplitude consistently with experimental ob-servations. A new prediction of our analysis is that, fora fixed field amplitude, the widths of the two plateausscale as Ω − and Ω − . We then numerically calculatethe HHC spectrum and show that multi-step plateausemerge when both the field amplitude and the staggeredpotential are strong enough.The rest of this paper is organized as follows. In Sec. II,we introduce our model Hamiltonian in the Floquet for-mulation. We also summarize the properties of the Flo- a r X i v : . [ c ond - m a t . o t h e r] N ov quet eigenstates and explain how to calculate the HHCspectrum. Section III summarizes the symmetry proper-ties of the HHC spectrum. In Sec. IV, we obtain an ana-lytic form of the HHC spectrum in the single-band limit,where the staggered potential is absent, and discuss theplateau in the spectrum. We also obtain the asymptot-ically exact Floquet eigenstates analytically for small ormoderate field amplitude. By using these analytic formsof eigenstates, we develop in Sec. V a perturbation the-ory with respect to the staggered potential, and discussthe new plateau induced by the potential. In Sec. VI, wenumerically analyze the HHC spectrum where the pertur-bation theory is not applicable. Section VII summarizesthe results with concluding remarks. In the Appendix,we provide supplemental technical details consolidatingthe discussions in the main text. II. FORMULATION OF THE PROBLEM
In this section, we derive formulas for the HHC interms of the Floquet eigenstates. We investigate a sit-uation in which a monochromatic ac electric field isturned on at time t = 0. We solve the time-dependentSchr¨odinger equation by invoking the Floquet eigenstatesand calculate the Fourier components of the electric cur-rent for the solution. A. Model
In this paper, we study the response of electron sys-tems on a lattice with unit cell containing multiple sites.As the simplest model, we consider a model of electronson a one-dimensional chain with the staggered potential,which doubles the size of the unit cell. We note that it isstraightforward to generalize the following arguments tothe cases of potentials with periodicity larger than twoand the results do not change qualitatively. The Hamil-tonian is given byˆ H = L (cid:88) j =1 (cid:104) t (ˆ c † j ˆ c j +1 + ˆ c † j +1 ˆ c j ) + Q ( − j ˆ c † j ˆ c j (cid:105) , (1)where ˆ c j (ˆ c † j ) is the annihilation (creation) operator forthe electron on the j -th site. We have ignored the spindegrees of freedom, which, if considered, only multiplythe following results by the factor of 2. The length of thechain is given by 2 L with an even number L ( > t denotes the transfer integral, Q is the amplitudeof the staggered potential, and we restrict ourselves tothe case of | Q | < Q splits the single cosineband into two. The physical systems described well by this model include some binary compounds with chemi-cal formula AB and the electrons in the presence of staticCDW order with period two. In the following, the unit ofenergy is fixed so that t = 1 /
2. This implies that the halfof the total band width for Q = 0 is set to unity in ourunit. Thus our unit of energy is read typically as 2.5 eVfor semiconductors and 0.25 eV for one-dimensional or-ganic conductors.We make a remark on the relationship between ourmodel (1) and the two-band models for typical semi-conductors. Generally speaking, multiple bands areformed from several single-electron states in the unit cell,and there are two typical cases. In the first case, theband multiplicity corresponds to the number of differentatomic orbitals at each site, and this is a standard setupfor semiconductors. In the second case, the band multi-plicity is the number of different sublattice sites in theunit cell, and we focus on this case in this paper. A tight-binding Hamiltonian can model both cases, and applyingthe electric field generally induces interband transitionsof electrons regardless of the origin of multiple bands, al-though their matrix elements depend on details such asthe type of atomic orbitals and the position of sublatticesites.The Hamiltonian (1) is diagonalized in the momen-tum space by the Fourier transformation for two-siteunit cells: ˆ a k = L − / (cid:80) Lj =1 e − i k (2 j ) ˆ c j and ˆ b k = L − / (cid:80) Lj =1 e − i k (2 j +1) ˆ c j +1 . Here the distance a be-tween the neighboring sites is set to unity, and the lat-tice momentum k takes the values of k = πn/L ( n = − L/ , − L/ , . . . , L/ − H = (cid:88) k ˆ φ † k H ( k ) ˆ φ k (2)with ˆ φ † k = (ˆ a † k ˆ b † k ) and the 2 × H ( k ) = cos k σ x + Qσ z , (3)where σ x and σ z are the Pauli matrices. The two eigen-values of H ( k ) are given by (cid:15) ± ( k ) = ± (cid:112) cos k + Q (4)and we refer to (cid:15) + ( k ) and (cid:15) − ( k ) as the upper and thelower bands, respectively. The band gap is given by2 | Q | , which is the energy difference at the Brillouin-zoneboundary k = − π/ E ( t )are taken into account in terms of the vector poten-tial A ( t ), which satisfies d A ( t ) / d t = − E ( t ). Through-out this paper, we assume that the electric field and,hence, the vector potential are homogeneous in space.The vector potential modifies the Hamiltonian (1)by the gauge-invariant Peierls substitution: ˆ c † j ˆ c j → ˆ c † j ˆ c j e i eA ( t )( j − j ) , where − e denotes the electron charge.Correspondingly, Eq. (2) is replaced byˆ H ( t ) = (cid:88) k ˆ φ † k H ( k, t ) ˆ φ k (5)where H ( k, t ) = H ( k + eA ( t )) = cos [ k + eA ( t )] σ x + Qσ z . (6)We note that this time-dependent Hamiltonian is diag-onal in k since the vector potential does not break thetranslation symmetry. B. Time evolution
In the present work, we consider the dynamics inducedby the ac electric field E ( t ) with frequency Ω ( > t = 0: E ( t ) = − E cos Ω t for t > E ( t <
0) = 0. This is represented by the followingvector potential: A ( t ) = A θ ( t ) sin Ω t ; A = E Ω . (7)The time dependence in the Hamiltonian (5) is now givenby H ( k, t ) = cos( k + F sin Ω t ) σ x + Qσ z (8)for t >
0. Here the dimensionless parameter F ≡ eA = eE Ω = Ω B Ω (9)quantifies the strength of the coupling to the input elec-tric field and Ω B = eE is the so-called Bloch frequency.Our monochromatic input (7) has two advantages.First, the high harmonics are well defined as multi-ples of Ω, in contrast to polychromatic inputs suchas a pulse [44]. Second, the time-dependent Hamilto-nian H ( k, t ) (8) becomes periodic in t > T = 2 π/ Ω. We will utilize this periodicity in thefollowing to solve the time-dependent Sch¨odinger equa-tion. We note that the input (7) has additional sym-metries A ( t ) = − A ( T / t ) = − A ( T − t ), which imply H ( − k, t ) = H ( k, t + T /
2) = H ( k, T − t ).As for the initial condition ( t < (cid:104) ˆ c † j ˆ c j (cid:105) + (cid:104) ˆ c † j +1 ˆ c j +1 (cid:105) = 1 in each unit cell, and the system is inthe ground state of ˆ H : | Ψ (cid:105) = (cid:89) k (cid:104) ˆ φ † k · (cid:126)ψ k (cid:105) | (cid:105) ; (cid:126)ψ k = (cid:18) ψ k,a ψ k,b (cid:19) , (10)where the product runs over all k ’s in the Brillouin zoneand | (cid:105) denotes the Fock vacuum [45]. Namely, (cid:126)ψ k cor-responds to the one-particle wave function with the neg-ative energy (cid:15) − ( k ). Then we are interested in the evolu-tion of the many-body state | Ψ( t ) (cid:105) = (cid:81) k [ ˆ φ † k · (cid:126)ψ k ( t )] | (cid:105) , in which each (cid:126)ψ k ( t ) obeys the one-particle Schr¨odingerequation i dd t (cid:126)ψ k ( t ) = H ( k, t ) (cid:126)ψ k ( t ) (11)with the initial condition (cid:126)ψ k ( t = 0) = (cid:126)ψ k . Here we haveused the fact that ˆ H ( t ) is diagonal in k and the particlenumber ˆ φ † k ˆ φ k for each k is conserved.Owing to the periodicity in time, the general solutionsof the time-dependent Sch¨odinger equation (11) are ob-tained by the Floquet theory [24]. The Floquet Hamil-tonian is given by [46] H F mn ( k ) = n Ω δ mn + (cid:90) T d tT H ( k, t )e − i( m − n )Ω t (12)= ( n Ω + Qσ z ) δ mn + i m − n J m − n ( F ) cos (cid:20) k − ( m − n ) π (cid:21) σ x (13)for each pair of integers m and n ( −∞ < m, n < ∞ ).Here is the 2 × J n ( z ) denotes theBessel function of the first kind. The Floquet eigenstates χ ( k ) = { (cid:126)χ n ( k ) } n are defined by (cid:88) n H F mn ( k ) (cid:126)χ αn ( k ) = E α ( k ) (cid:126)χ αm ( k ) , (14)and E α ( k )’s are the Floquet eigenvalues. The corre-sponding time-dependent physical state is given by (cid:126)χ α ( k, t ) = e − i E α ( k ) t (cid:88) n (cid:126)χ αn ( k )e i n Ω t . (15)This becomes a solution of Eq. (11) and satisfies (cid:126)χ α ( k, t + T ) = e − i E α ( k ) T (cid:126)χ α ( k, t ).Although the Floquet Hamiltonian (12) has an infi-nite number of eigenstates, only two of them are phys-ically independent and this number is the dimensionof H ( k, t ). This is because, if { (cid:126)χ n ( k ) } n is a Floqueteigenstate with eigenvalue E ( k ), then, shifting this inthe Floquet space by any integer M leads to anothereigenstate { (cid:126)χ n + M ( k ) } n with eigenvalue E ( k ) − M Ω, andthese shifted eigenstates all describe the same evolu-tion (15). To avoid this redundancy, we may choose thetwo Floquet eigenstates with eigenvalues in the interval[ − Ω / , Ω / χ α ( k ) = { (cid:126)χ αn ( k ) } n [ α = 1 , E ( k ) ≤ E ( k )] be the inequivalent Floquet eigen-states thus obtained, which are normalized and orthogo-nal to each other χ α ( k ) † χ β ( k ) = δ αβ .Equation (11) with our initial condition is solved byexpanding the initial state (cid:126)ψ k in terms of (cid:126)χ α ( k, t = 0) = (cid:80) n (cid:126)χ αn ( k ) ≡ (cid:126)X α ( k ). This expansion is always possiblesince (cid:80) α (cid:126)X α ( k ) (cid:126)X α † ( k ) = , and the expansion coeffi-cients are calculated as w α ( k ) = (cid:126)X α ( k ) † (cid:126)ψ k , which satisfy -1-0.5 0 0.5 1 0 0.25 0.5 -0.125 0 0.125 0 0.25 0.5-1-0.5 0 0.5 1 0 0.25 0.5-0.5 0 0.5 0 0.25 0.5 (a)(b)(c) (d)(e)( f ) F l oque t E i gene v a l ue F l oque t E i genen v a l ue -1-0.5 0 0.5 1 0 0.25 0.5 -1-0.5 0 0.5 1 0 0.25 0.5 FIG. 1. (a) Floquet eigenvalues for F = 0 and Ω = 1 .
0. Thesolid lines correspond to the energy bands (cid:15) ± ( k ) [Eq. (4)], andthe dashed lines to the other Floquet bands. The shaded arearepresents the representative interval [ − Ω / , Ω / F = 0 .
5, (b)the representative Floquet eigenvalues and (c) the weight dif-ference | w ( k ) | − | w ( k ) | are plotted against k . The panels(d)-(f) are similar plots to (a)-(c), where Ω is changed to 0 . Q = 0 . (cid:80) α | w α ( k ) | = 1. Then, the solution of Eq. (11) is givenby (cid:126)ψ k ( t ) = (cid:88) α w α ( k )e − i E α ( k ) t (cid:88) n (cid:126)χ αn ( k )e i n Ω t . (16)We note that Eq. (16) holds true only for t ≥
0, and (cid:126)ψ k ( t <
0) = e − i (cid:15) − ( k ) t (cid:126)ψ k for t ≤ H F mn ( k ) with asufficiently large cutoff for | m | and | n | , and the expan-sion coefficients w α ( k ) are calculated from them. Fig-ure 1 shows in the left column the results for Ω = 1.First, the panel (a) shows the Floquet eigenvalues for F = 0 for reference, i.e. (cid:15) ± ( k ) − M Ω for some M ’s. Thepanel (b) shows the representative Floquet eigenvalues E α ( k ) for F = 1. They change smoothly with F atmost k ’s, but the crossing points of the different bandsat F = 0 become anticrossings at F (cid:54) = 0. The panel (c)shows the weight difference | w ( k ) | − | w ( k ) | of the twoFloquet eigenstates. It changes sign at the anticrossingpoints, and its modulus reduces also at the Brillouin-zoneboundary. Figures 1(d)-(f) show the corresponding datafor Ω = 0 .
25. The lower input frequency increases thenumber of anticrossing points in the Floquet bands, andthis results in more oscillations in the weight difference.
C. High-harmonic current (HHC)
Here we use the solution (16) represented by the Flo-quet eigenstates (14) to calculate the time evolution ofthe electric current, which is the source of radiation. Wewill show that the current spectrum consists of a discretepart peaked at n Ω ( n ∈ Z ) and a continuous part. Theformer, the high-harmonic current, works as the sourceof HHG.The electric current density is obtained as the expec-tation value of the operatorˆ I ( t ) ≡ L ∂ ˆ H ( t ) ∂A ( t ) = 12 L (cid:88) k ˆ φ † k I k ( t ) ˆ φ k , (17)which is again diagonal in k as seen from Eq. (5). Since A ( t ) is periodic in time, I k ( t ) is also periodic and, hence,expanded in a Fourier series as I k ( t ) = ∂H ( k, t ) ∂A ( t ) = − e sin( k + F sin Ω t ) σ x = (cid:88) n I k,n e − i n Ω t . (18)We note that I k, − n = I † k,n since I k ( t ) is Hermitian. The2 × I k,n is given by I k,n = − i n eJ n ( F ) sin (cid:16) k + nπ (cid:17) σ x ≡ − ev n ( k ) σ x . (19)Equation (19) follows from Eqs. (5) and (3) and the her-miticity implies I k, − n = ( − n I k,n = I ∗ k,n . Evaluatingthe expectation value of Eq. (17) for the solution (16),we obtain I ( t ) = (cid:104) Ψ( t ) | ˆ I ( t ) | Ψ( t ) (cid:105) = 12 L (cid:88) k (cid:126)ψ k ( t ) † I k ( t ) (cid:126)ψ k ( k ) (20)= 12 L (cid:88) k,α,β w α ( k ) ∗ w β ( k ) (cid:88) m,n,l (cid:126)χ αm ( k ) † I k,l (cid:126)χ βn ( k ) × e − i[ E β ( k ) − E α ( k )+( m + l − n )Ω] t . (21)Equation (20) consists of two kinds of contributions,which are diagonal ( α = β ) and off-diagonal ( α (cid:54) = β ) interms of the labels for the Floquet eigenstates [47]: I ( t ) = I H ( t ) + I C ( t ) , (22) I H ( t ) ≡ (cid:88) N I H ( N )e − i N Ω t , (23) I C ( t ) ≡ (cid:90) ∞−∞ d ω π I C ( ω )e − i ωt . (24)The part I H ( t ) corresponds to the diagonal part, and wehave I H ( N ) ≡ L (cid:88) k,α | w α ( k ) | I α H ( k, N ) , (25) I α H ( k, N ) ≡ (cid:88) n,l (cid:126)χ αn − l + N ( k ) † I k,l (cid:126)χ αn ( k ) . (26)This contribution only involves a discrete set of the har-monics of the input frequency Ω. On the other hand, I C ( t ) involves the k -dependent frequencies E α ( k ) − E β ( k ),which lead to a continuous spectrum I C ( ω ) in the ther-modynamic limit, L → ∞ . In the time domain, I C ( t )decays [48] whereas I H ( t ) remains to oscillate as t → ∞ ,and the harmonic part dominates for sufficiently large t .Thus we regard Eq. (25) as the source of the N -thorder HHG and refer to I H ( N ) as the HHC spectrum.According to the classical electromagnetism, the totalradiation power P rad ( N ) from the N -th HHC is propor-tional to ( N Ω) | I H ( N ) | . In the following, we investigatethe symmetry aspects and the N -dependence of I H ( N )rather than P rad ( N ) for comparison to the related the-oretical studies. One can immediately obtain P rad ( N )from I H ( N ) if necessary, and the plateaus discussed be-low will become clearer when plotted for P rad ( N ) due tothe factor ( N Ω) .We note that Eq. (25) is a generalization of similarformulas in the literature (see e.g. Eq. (28) in Ref. [36]).In the literature, it is assume that some of the Floquetbands are fully occupied and the quantity I α H (cid:48) ( N ) = 12 L (cid:88) k I α H ( k, N ) (27)is discussed. On the other hand, our formula (25) in-volves the weight factor | w α ( k ) | on each Floquet eigen-state, and | w α ( k ) | can take any value between 0 and1. The effects of the fractional weight factor becomemore significant for a stronger electric field or near theBrillouin-zone boundary and the anticrossing points asshown in Fig. 1. III. SYMMETRY ASPECTS
Formulas (25) and (26) for the HHC based on theFloquet eigenstates have the advantage that the conse-quences of the symmetry are manifest. In this section,we first reproduce the important known property that I α H (cid:48) ( N ) (27) vanishes for any even N owing to the in-version symmetry (see e.g. Ref. [29]). Then we discussthe symmetry of the weights | w α ( ± k ) | and show that I H ( N ) (25) vanishes for any even N in our choice of thevector potential A ( t ), although this does not hold oncethe initial phase of the input filed is shifted.We begin by noting that the inversion symmetry H ( − k ) = H ( k ) (28)breaks down at time t > H ( − k, t + T /
2) = H ( k, t ) owing to A ( t + T /
2) = − A ( t ). This leads to thefollowing symmetry for the Floquet Hamiltonian H F mn ( − k ) = ( − m − n H F mn ( k ) , (29) which implies (cid:126)χ αn ( − k ) = ( − n (cid:126)χ αn ( k ) for an appropri-ate choice of the overall phase. Together with I − k,n =( − n +1 I k,n , it follows from Eq. (26) I α H ( − k, N ) = ( − N +1 I α H ( k, N ) . (30)This means I α H (cid:48) ( N ) = 0 for even N since the contribu-tions from ± k cancel out with each other. We remarkthat Eq. (26) holds true also for any inputs as long as A ( t + T /
2) = − A ( t ) is satisfied.The inversion symmetry between the weights | w α ( ± k ) | follows from yet another symmetry H ( − k, T − t ) = H ( k, T + t ), which leads to H F mn ( − k ) = H F nm ( k ) = H F mn ( k ) ∗ . (31)This implies (cid:126)χ αn ( − k ) = (cid:126)χ αn ( k ) ∗ and (cid:126)X α ( − k ) = (cid:126)X α ( k ) ∗ with appropriate choices of the overall phases. Besides,since Eqs. (28) and (3) ensure (cid:126)ψ − k = (cid:126)ψ k = (cid:126)ψ k ∗ , weobtain w α ( − k ) = w α ( k ) ∗ and, hence, | w α ( − k ) | = | w α ( k ) | (32)for all α ’s. The inversion symmetry of | w α ( ± k ) | (32), to-gether with Eq. (30), means that the HHC spectrum (25)vanishes for even N ’s.We note that this property (32) is violated if we shiftthe initial phase of the input as A ( t ) = A sin(Ω t + θ ).When 0 < θ < π , the time evolution for ± k occurasymmetrically, and we have (cid:126)χ αn ( − k ) = (cid:126)χ αn ( k ) ∗ e nθ and | w α ( − k ) | (cid:54) = | w α ( k ) | . This asymmetry in the weightsof Floquet eigenstates leads to I H ( N ) (cid:54) = 0 for even N ’seven if ˆ H has the inversion symmetry. In the following,we restrict ourselves to the case of θ = 0 and discuss the N -dependence of the HHC spectrum for odd N ’s. IV. SINGLE-BAND LIMIT
In this section, we discuss the case where Q = 0 and thestaggered potential is absent. We refer to this case as thesingle-band limit because we can also use the single-siteunit cell and the Brillouin zone is doubled, where each k has only one energy band cos k . However, to comparewith the case of Q (cid:54) = 0, we keep using the two-site unitcell. As discussed in Ref. [49], the HHC are present due tothe nonlinearity of the Peierls substitution, in contrast tothe continuous models (see e.g., Ref. [33]). In Sec. IV A,we calculate the HHC spectrum for Q = 0 in our for-mulation and obtain results consistent with Ref. [49]. InSec. IV B, we derive the asymptotically exact eigenstatesfor the Floquet Hamiltonian for small | F | ( (cid:46) Q (cid:54) = 0 case inSec. V. A. The HHC spectrum
In this special case of Q = 0, the time-evolution oper-ator for Eq. (11) is exactly obtained since the Hamilto-nians at different times commute with each other. In the2 × U k ( t ) = exp (cid:20) − i σ x (cid:90) t d τ cos( k + F sin Ω τ ) (cid:21) . (33)It is noteworthy that this commutes with the cur-rent matrix [Eqs. (18) and (19)], and therefore U k ( t ) † I k ( t ) U k ( t ) = I k ( t ).It is straightforward to calculate the HHC spectrumfrom Eqs. (18) and (19). Noting that the initial state (cid:126)ψ k is the eigenvector of σ x with − I H ( k, N ) = e i N J N ( F ) sin (cid:18) k + N π (cid:19) (34)Its k -sum vanishes for even N ’s, and this is consistentwith the inversion symmetry as discussed in Sec. III. Forodd N ’s, we obtain I H ( N ) = lim L →∞ L (cid:88) k I H ( k, N ) = i eπ J N ( F ) . (35)We note that, for an initial state at arbitrary filling, thesum over k in Eq. (35) is restricted, and the result ismultiplied by sin( πρ ), where ρ is the electron density and1 / I ( t ) = eπ sin ( F sin Ω t ) (36)in the limit of L → ∞ . One can check that the Fourierexpansion of Eq. (36) reproduces Eq. (35). Furthermore,Eq. (36) implies that the HHC spectrum contains onlyharmonics of Ω and the continuous part does not exist, I C ( t ) = 0, in the single-band limit.Equation (35) implies that the HHC spectrum qualita-tively changes depending on whether | F | < | F | > | J N ( F ) | with fixed F isapproximately constant for | N | (cid:46) | F | and rapidly decaysfor | N | (cid:38) | F | . Thus, when | F | <
1, the HHC spectrummerely decays as | N | increases. On the other hand, when | F | >
1, a plateau emerges in the spectrum and its widthis given by | F | , which is proportional to E and Ω − .These features are shown in Fig. 2. This plateau is es-sentially the same as the one discovered by Pronin andcoworkers [49].We note, however, that this plateau is too narrow toexplain the experiment by Ghimire et al. [1] that detectedthe harmonics up to the 25th at F ∼
5. They also showedthe presence of a wider plateau if the band dispersion isdeformed from cos k even in the single-band case. Laterin Sec. V, we will show another mechanism for a widerplateau, i.e., the staggered potential Q splitting a singleband into two. This wider plateau has a different scalingof its width with input frequency Ω. -8 -4
0 5 10 15 20
Harmonic Order
FIG. 2. The HHC spectrum in the single-band limit Q = 0[Eq. (35)]. Each data set corresponds to the field strength F =0 . B. Floquet Eigenstates at | F | < We will show in Sec. V that the staggered potential Q produces another plateau even when | F | <
1. For thispurpose, we here derive the Floquet eigenstates χ α ( k ) inthe single-band limit.To obtain the Floquet eigenstates, we first calculate theFourier expansion of the time-evolution operator U k ( t ).Since the exact result is very complicated, we considerthe case of | F | < k + F sin Ω t ) = (cid:80) n J n ( F )Re[e i( k + n Ω t ) ] in Eq. (33)by its partial sum of − ≤ n ≤
1. Within this approxi-mation, the time-evolution operator is given by U T k ( t ) = e − i tJ ( F ) cos kσ x e − i z k [cos(Ω t ) − σ x (37)= e − i[ tJ ( F ) cos k − z k ] σ x (cid:88) n ( − i σ x ) n J n ( z k )e i n Ω t (38)with z k ≡ t J ( F )Ω sin k ∼ t F Ω sin k. (39)Here the superscript T indicates that the truncation isperformed, and we have recovered the transfer integral t , which have been set to 1 / J n ( F ) with | n | ≥
2, which are O ( F ), areneglected in this approximation.The solutions of the time-dependent Schr¨odinger equa-tion is immediately obtained within this approximationsince U T k ( t ) contains only σ x . By applying U T k ( t ) ontothe eigenstates of σ x , σ x (cid:126)ζ ± = ± (cid:126)ζ ± i.e. (cid:126)ζ ± = (cid:0) ± (cid:1) / √ U T k ( t ) (cid:126)ζ α = e i αz k − i E α ( k ) t (cid:88) n ( − i α ) n J n ( z k )e i n Ω t (cid:126)ζ α (40)for α = ± , where we have introduced E α ( k ) = αJ ( F ) cos k. (41)Comparing Eqs. (15) and (40), one finds that the twoFloquet eigenvalues are E ± ( k ) and their eigenstates aregiven by (cid:126)χ ± n ( k ) = ( ± i) n √ J n ( z k ) (cid:18) ± (cid:19) ≡ c ± n ( k ) (cid:126)ζ ± . (42)Here we have ignored the phase factor e ± i z k and thiscorresponds to the choice of the global phase of (cid:126)χ α ( k, t ).We note that we do not require E α ( k ) ∈ [ − Ω / , Ω /
2) inthe analytical calculations for convenience.The above approximation is equivalent to truncatingthe off-diagonal elements in the Floquet Hamiltonian as H T,F mn ( k ) ≡ (cid:40) H F mn ( k ) ( | m − n | ≤ . (43)The complete set of the eigenstates χ α,M ( k ) of the trun-cated Floquet Hamiltonian H T,F ( k ) are defined for inte-ger M ’s as (cid:126)χ α,Mn ( k ) = (cid:126)χ αn + M ( k ) = c αn + M ( k ) (cid:126)ζ α . (44)Then they satisfy the eigenvalue equation (cid:88) n H T,F mn ( k ) (cid:126)χ α,Mn ( k ) = [ E α ( k ) + M Ω] (cid:126)χ α,Mm ( k ) . (45)Now we discuss the distribution of a Floquet eigenstateover the Floquet space, or index n : p αn ( k ) ≡ | (cid:126)χ αn ( k ) | = J n ( z k ) . (46)As we have seen above, this distribution is approximatelyconstant for | n | (cid:46) | z k | and rapidly decays for | n | (cid:38) | z k | .Therefore, if Ω is smaller enough than t , a plateau withwidth | z k | emerges in p αn ( k ) even for F < p αn ( k ) is larger than that in the HHC bythe factor ξ = | z k || F | ∼ t Ω , (47)for Ω < t , where | z k | is the k -space average of | z k | .We note that, in the single-band limit, the plateau in p αn ( k ) has nothing to do with the HHC spectrum. Thisis because the time-dependent state (40) is always pro-portional to either of (cid:126)ζ α ’s and the high-harmonic termsin Eq. (40) amount to an overall phase factor. In fact,we have shown in Sec. IV A that the HHC spectrum doesnot show a plateau for | F | <
1, although p αn ( k ) can showa plateau. We will show in Sec. V, however, that theplateau in p αn ( k ) is converted into the HHC spectrumonce the staggered potential is turned on.We remark on the work in Ref. [38] that studied thecase of quasistatic input field with frequency Ω muchsmaller than the Bloch frequency Ω B . This correspondsto the limit of F = Ω B / Ω → ∞ in the present study, and the Bloch oscillation occurs many times within one periodof input time dependence. Although this differs fromthe typical situation in the present study, one can applythe present formulation without any problem also to thisparameter regime, as far as the input field is periodicin time. The harmonics in output are multiples of Ω,which distribute densely in the frequency space for smallΩ and one of them is very close to the Bloch frequency, N B Ω ∼ Ω B . Considering that the Floquet Hamiltonianhas the largest matrix element for this harmonics N B ,we expect that the output spectrum shows peaks aroundmultiples of N B Ω, which is consistent with the result ofRef. [38] predicting harmonics of Ω B . For this regime,one needs to diagonalize the Floquet Hamiltonian witha very large dimension greater than N B , and we do notfurther analyze this case. V. NEW PLATEAU INDUCED BYSTAGGERED POTENTIAL
In this section, we study the case of Q (cid:54) = 0 and ex-amine the effects of the staggered potential on the HHCspectrum. We will develop an analytical perturbativeapproach to the effects of the staggered potential Q , andmainly focus on the region of | F | < Q = 0 limit does notshow a plateau, we will show that a plateau appears inthe HHC spectrum at the order of Q .Since | F | (cid:46)
1, we may use the truncated FloquetHamiltonian H T,F ( k ) (43). We expand its eigenstatesas a polynomial in Q : χ α,M ( k ) + Q ν α,M ( k ) + Q λ α,M ( k ) + · · · , (48)where χ α,M ( k ) is the Floquet eigenstates in the single-band limit ( Q = 0) discussed in Sec. IV B. Correspond-ingly, we expand the HHC spectrum for each Floqueteigenstate (26) also as a polynomial in Q : I α, ( k, N ) + QI α, ( k, N ) + Q I α, ( k, N ) + · · · , (49)and we will calculate I α, ( k, N ) and I α, ( k, N ) below. A. Perturbation Theory
The perturbation to the Floquet Hamiltonian is V F mn = Qσ z δ mn (50)and this interchanges the two eigenstates (cid:126)ζ ± of σ x . Itsmatrix elements between the unperturbed eigenstates aregiven by χ α,M ( k ) † V F χ β,M (cid:48) ( k ) = Q ( α i) M − M (cid:48) J M − M (cid:48) (2 z k ) δ α, − β . (51)Let us focus on the correction for M = 0 owing to thephysical equivalence of the Floquet eigenstates. The first-order correction ν α, ( k ) is obtained by the standard first-order perturbation theory as ν α, ( k ) = (cid:88) M (cid:54) =0 B αM ( k ) χ − α,M ( k ) (52)with B αM ( k ) = ( − α i) M J M (2 z k ) M Ω + 2 αJ ( F ) cos k . (53)We note that |B αM ( k ) | also shows a plateau structure in M due to the Bessel function. When the vanishing of thedenominator of Eq. (52) is ignored, the M -dependence of |B αM ( k ) | is approximately constant for | M | (cid:46) | z k | andrapidly decays for larger | M | as shown in Fig. 3(a).Although the correction ν α, ( k ) spreads over variousFloquet bands, the HHC does not change at the firstorder of Q , or the first-order contribution to the HHCvanishes: I α, ( k, N ) = 0 . (54)This is because the current matrix (19) is proportionalto σ x and χ α, † σ x ν α,M = 0.We make a remark on the vanishing of the denom-inator in Eq. (52) for some k . This condition impliesthe resonance between two Floquet eigenstates and oneneeds a degenerate perturbation analysis. We show inAppendix A that this resonance actually gives contri-butions of O ( Q ). However, this contribution has thesame N dependence as the HHC spectrum in the single-band limit (34), and, hence, does not show a plateau for | F | (cid:46) O ( Q ) correction of the HHC spec-trum due to the matrix elements between ν α,M ’s. Thisis a part of I α, ( k, N ) and we define this as I α, A H ( k, N ) ≡ (cid:88) n,l (cid:126)ν α, n − l + N ( k ) I k,l (cid:126)ν α, n ( k ) . (55)In fact, another contribution of O ( Q ) comes from thesecond-order correction of the wave function λ α,M ( k ). InAppendix B, we show that its N -dependence is similar tothat of Eq. (55). By invoking Eqs. (52) and performingsome algebra, we obtain I α, A H ( k, N )= αe (cid:88) M,M (cid:48) B αM ( k ) ∗ B αM (cid:48) ( k ) v M − M (cid:48) + N ( k ) . (56)Equation (56) implies that, when 4 | z k | >
1, the HHCspectrum I α, A H ( k, N ) shows a plateau for | N | (cid:46) | z k | andrapidly decays for larger | N | as understood as follows.Since we are considering | F | (cid:46) v n ( k ) ∝ J n ( F )rapidly decays as | n | increases, the sum over M (cid:48) is dom-inated by M (cid:48) = N + M and Eq. (56) is approximated -15 -11 -7 -3
0 10 20 30 all
Harmonic Order (a) (c) -6 -3 -20 -10 0 10 20 (b) FIG. 3. (a) Typical behavior of |B αM ( k ) | [Eq. (53)] calculatedfor F = 0 .
5, Ω = 0 . α = +, and k = 7 π/
20. The plateauregion M (cid:46) z k is shown by arrows. (b) Schematic illus-tration for the overlap between B αM ( k ) and B αN + M ( k ). For N (cid:38) z k , the overlap rapidly decays. (c) The second-ordercorrection I α, ( k, N ) to the HHC spectrum (55) for k = 7 π/ k = ( π/ m/
10) with m = 0 , , . . . ,
9. The arrow indicates 4 z k for k = 7 π/ as αev ( k ) (cid:80) M B αM ( k ) ∗ B αN + M ( k ). As Fig. 3(b) shows,this sum, or the overlap between B αM ( k ) and B αN + M ( k ),rapidly decays for | N | (cid:38) | z k | , whereas it changes ratherslowly for | N | (cid:46) | z k | [50]. This is how a plateau appearsin the HHC spectrum at the second order of the staggeredpotential Q . The above argument is confirmed by the nu-merical results shown in Fig. 3(c), where I α, A H ( k, N ) iscalculated as in Eq. (55) and plotted for a representative k and the average over k .As a result of these analyses, we propose N cut = 8 π | J ( F ) | t Ω (57)as an indicator of the plateau width, or the high-energycutoff order of the HHC spectrum. Here we have recov-ered the transfer integral t , which has been set to 1 / (cid:80) k I α, A H ( k, N ) and4 z k depends on k . Averaging | z k | over the Brillouin zone,we obtain Eq. (57). Since the averaging is a crude ap-proximation, the numerical factor 8 /π in Eq. (57) shouldnot be taken very seriously.Equation (57) can be used to derive the onset fieldstrength F onset at which the plateau sets in. Whether aplateau exists or not should correspond to N cut (cid:46) N cut (cid:38)
1, respectively. Thus F onset is estimated by thecondition N cut = 1, which leads to F onset = Ω t , (58)where we have ignored the numerical factor π/ J ( F ) by F/ F onset is smallenough. We emphasize that F onset can be less than 1if Ω < t . B. Scaling properties of the plateau
Now we discuss how N cut depends on the amplitude E and the frequency Ω of the input ac electric field. When | F | (cid:46)
1, we approximate J ( F ) (cid:39) F/ N cut ∼ | eaE t | ( (cid:126) Ω) , (59)where we have used Eq. (9) and recovered a and (cid:126) thathave been set to unity.Equation (59) shows that the cutoff order N cut is pro-portional to the amplitude E rather than the power E of the input electric field. This is consistent with the ex-perimental observations [1] and a unique feature of theHHG in solids in contrast to that in gases.A remarkable prediction of Eq. (59) is that the cutofforder N cut is proportional to Ω − rather than Ω − for afixed field amplitude | E | . This originates from the in-trinsic property of the Floquet eigenstate χ α,M ( k ). Asshown in Sec. IV B, this eigenstate distributes over theFloquet index and the width of the distribution is propor-tional to | z k | ∼ J ( F ) / Ω ∼ F/ Ω ∝ Ω − since F = Ω B / Ω.Thus the scaling N cut ∝ Ω − is a signature of the Floqueteigenstate, which could be tested in experiments.We remark that the cutoff energy defined by E cut ≡ N cut (cid:126) Ω ∼ | eaE t | (cid:126) Ω (60)has a slightly different scaling. The cutoff energy is pro-portional to E and Ω − . This scaling could also betested experimentally if several frequencies for the inputare available. C. Numerical Verification
Let us numerically verify the above analytical argu-ments based on perturbation theory. Figure 4 shows theHHC spectrum I H ( N ) calculated as in Eq. (25) for sev-eral parameter sets ( F, Ω) with L = 10 and Q = 0 . H F mn ( k ), and the cutoff for the Floquet index hasbeen chosen as − ≤ m, n ≤ F, Ω) inFig. 4. For the F = 1 . N cut as 13 (Ω = 0 . . .
1) as indicated by arrows in the figure. Theratios between these numbers are in good agreement withthe analytical prediction (57), which states that N cut isproportional to Ω − with fixed F . We also assign N cut =11 for the ( F, Ω) = (0 . , .
1) data, which is approximatelyconsistent with Eq. (57).We have also verified that the Q -dependence of themagnitude | I H ( N ) | is consistent with the perturbationanalysis. This requires a careful treatment due to theresonances between the Floquet eigenstates as remarked -15 -12 -9 -6 -3
0 5 10 15 20 25 30
Harmonic Order
FIG. 4. The HHC spectrum calculated for the four parametersets ( F, Ω) as indicated in the figure. The parameter Q isset 0 .
01. The arrow indicates the end of the plateau that isidentified from the plot. in Sec. V A. We show the details of the verification inAppendix A.
VI. BEYOND PERTURBATION
In Secs. IV and V, we have analytically investigated theHHC spectrum in the two limiting cases: (i) Q = 0 andarbitrary F , and (ii) | F | (cid:46) | Q | . In the othercases, we numerically calculate the HHC spectrum andshow that the scalings in the limiting cases still hold ifeither | F | or | Q | is small, whereas a qualitative differencesets in when both | F | and | Q | become large.First, we investigate the case of | F | ≥ | Q | . Figure 5(a) shows the HHC spectrum for F = 1 . . Q = 0 .
01. For F = 5 . ≤ N ≤ N ≥ F . The second plateau already existsat F = 1 .
0, and it is induced by the staggered poten-tial. We note that the width of the second plateau doesnot necessarily follow Eq. (57) since the truncation of theFloquet Hamiltonian is no longer justified for F (cid:38) | F | (cid:46) | Q | .Figure 5(b) shows the HHC spectrums for Q = 0 . . F, Ω) = (1 . , . Q = 0 .
01, for which the perturbation analysisworks well. For Q = 0 . .
2, the perturbation theoryis no longer justified since Q is not the smallest parame-ter, but its results are still valid approximately. Namely,compared with the data for Q = 0 .
01, the plateau hasalmost the same width, and its height is enhanced abouttwo orders of magnitude. Thus we conclude that thepresence of the plateau induced by the staggered poten-tial is not restricted to the region | Q | < Ω, but can be0
Harmonic Order Harmonic Order Harmonic Order -9 -6 -3
0 10 20 30 40 (a) (b) (c) -9 -6 -3
0 5 10 15 20 25 30 10 -15 -10 -5
0 10 20 30 40 50 60
FIG. 5. Numerically calculated HHC spectrum in the nonperturbative regimes including the three cases (a) small Q but large F , (b) small F but large Q , and (c) large F and Q . In all panels, Ω = 0 .
1. In panel (a), the three data sets correspond todifferent field strengths F = 5 . . . Q = 0 .
01. In panel (b), the staggered potentialstrength is varied: Q = 0 .
20 (circle), 0.10 (square), and 0.01 (triangle) with F = 1 .
0. In panel (c), the parameter set is (
F, Q ) =(3.0,0.30) (filled circle), (1.0,0.30) (open circle), (3.0,0.01) (filled square), and (1.0,0.01) (open square). Dashed lines indicatethe multi-step plateaus. extended to | Q | > Ω for | F | (cid:46) | F | nor | Q | is small. Figure 5(c) shows that the HHC spectrum for Q = 0 . F = 3 . N (cid:38) F = 1 .
0. Moreover, in addi-tion to the plateau in N (cid:46)
20, multi-step plateaus emergeas indicated by dashed lines in the figure. We also plotthe data for Q = 0 .
01 in Fig. 5(c) for comparison. Al-though the enhancement of the magnitude with F is com-mon for both Q ’s, the multi-step plateaus appear only forthe larger Q . Thus a qualitative difference arises in theHHC spectrum when neither | F | nor | Q | is small, andan approach beyond perturbation is desired to reveal theorigin of the multi-step plateaus.We make a remark on the possible relevance of theplateaus in N (cid:38)
30 in experiments. Ndabashimiye andcoworkers [5] have observed the HHG in rare-gas solidsand reported that a new plateau emerges around N ∼ VII. CONCLUSIONS
We have considered the setup where an ac electric fieldwith frequency Ω is turned on at some time. In this setup,the harmonics are well defined as multiples of Ω, and theHHC spectrum have been related to the Floquet eigen-states [see Eq. (25) and (26)]. Our formulation is a gener-alization of similar formulas in the literature because thecondition that the initial state is the ground state is takeninto account by the weights of each Floquet eigenstate. Inthis formulation, analytical approaches are feasible, andthe consequences of symmetries of the Hamiltonian areeasily tractable from the Floquet eigenstates and theirweights.On the basis of this formulation, we have investigated the HHC spectrum of electrons on a one-dimensionalchain with the staggered potential Q , which splits thesingle cosine band into two. In the single-band limit( Q = 0), we have confirmed the result [49] that a plateauof width | F | appears owing to the nonlinearity of thePeierls substitution for strong field | F | >
1, where the di-mensionless parameter F = Ω B / Ω ∝ E / Ω quantifies thecoupling between the electron and the electric field. Ournew finding is that the staggered potential Q induces an-other wider plateau emerging from weaker field | F | < Q = 0, we have shown that thewidth of plateau induced by the staggered potential is ξ | F | = | eaE t | / ( (cid:126) Ω) , which is proportional to E andlarger by the factor ξ = t / Ω than that in the absence ofthe staggered potential. Our result also provides a newprediction that the width of the plateau scales as Ω − .Since this scaling originates from the Floquet eigenstates,it could also be an experimental signature in identifyingthose states.We have numerically confirmed that our analytical re-sults hold qualitatively as far as either the field amplitude | F | or the staggered potential | Q | is small. This condi-tion includes the case where Ω is smaller than the bandgap. We have also analyzed the case where both | F | and | Q | are large and found that a more complex structurein the HHC spectrum involving the multi-step plateaus,which might be relevant in interpreting experiments andmerits further systematic studies.We make a remark on an implication on the HHG ex-periments in CDW materials. The cases Q = 0 and Q (cid:54) = 0 correspond to the phases above and below thetransition temperature T c of the CDW order. Because ofhigh carrier density above T c , the largest | E | is limitedin experiments to avoid damaging the samples and only afew harmonics would be observable. In fact, the 9th har-monic has been the highest observed in 2H-NbSe above T c [51]. Our results imply that, even with the same lim-ited | E | , the HHG is enhanced below T c and a plateau1could be observable in the spectrum. In this situation,our results also predict that the cutoff order N cut scalesas Ω − and Ω − above and below T c , respectively.As a concluding remark, we should mention that sev-eral physical processes are not taken into account in ourformulation. First, the correlation effects [23, 52–54] andthe order parameters [55] have been neglected in ourmodel. Second, the energy dissipation to the phononthermal bath [56]. has not been considered either, whichmight be relevant if the driving frequency is as low as1THz ∼ − . For the carrier bath, the Keldysh for-malism has been employed to calculate the HHC [57–60]. Third, we have treated the electric field classically.Quantum processes lead to spontaneous emission of pho-tons [61], which is not discussed in the present work.These effects are all important and our simple modelcould serve as a starting point in interpreting the ex-perimental data. ACKNOWLEDGEMENTS
Fruitful discussions with Y. Kayanuma, K. Shimo-mura, T. Tamaya, and K. Uchida are gratefully acknowl-edged. This work was supported by JSPS KAKENHI Grants No. JP16H06718 and JP18K13495. K.C. ac-knowledges financial support provided by the AdvancedLeading Graduate Course for Photon Science at the Uni-versity of Tokyo.
Appendix A: Resonances of Floquet eigenstates
In this appendix, we refine the perturbation theory inSec. V A by taking account of the resonances of the Flo-quet eigenstates, and calculate the O ( Q ) contributions tothe HHC spectrum.Let us take a positive integer (cid:96) satisfying (cid:96) < | J ( F ) / Ω | . Then there exist a momentum k (cid:96) such that E + ( k (cid:96) ) − (cid:96) Ω = E − ( k (cid:96) ) . (A1)This means that when Q = 0 the Floquet eigenstates χ + ,(cid:96) ( k (cid:96) ) and χ − , ( k (cid:96) ) are degenerate. We investigatethe Floquet eigenstates for Q (cid:54) = 0 in the vicinity of k = k (cid:96) , where the mixing of these states cannot be treatedby perturbation theory with respect to Q , and must betreated exactly.Since we are interested in the vicinity of k = k (cid:96) , the2 × χ + ,(cid:96) ( k (cid:96) ) and χ − , ( k (cid:96) ) is linearized in ∆ k ≡ k − k (cid:96) ,and we obtain H F,res ( k ) = (cid:18) E + ( k ) − (cid:96) Ω Q i (cid:96) J (cid:96) (2 z k ) Q ( − i) (cid:96) J (cid:96) (2 z k ) E − ( k ) (cid:19) ∼ E − ( k (cid:96) ) + (cid:18) − a i (cid:96) b ( − i) (cid:96) b a (cid:19) , (A2)where is the unit matrix, z k = 2 J ( F ) sin k/ Ω,and we have defined a ≡ u (cid:96) ∆ k , b ≡ ( c (cid:96) + d (cid:96) ∆ k ) Q , u (cid:96) ≡ J ( F ) sin k (cid:96) , c (cid:96) ≡ J (cid:96) (2 z k (cid:96) ), and d (cid:96) ≡ J (cid:48) (cid:96) (2 z k (cid:96) ) (cid:96)J ( F ) /J ( F ).The two eigenvalues of the linearized Hamiltonian (A2)are given by E − ( k (cid:96) ) ± ∆ E (∆ k ) (A3)with the energy splitting∆ E (∆ k ) ≡ (cid:113) Q ( c (cid:96) + d (cid:96) ∆ k ) + u (cid:96) (∆ k ) , (A4) and the two-fold degeneracy is lifted by the coupling Q .We emphasize that the energy splitting is not minimal at∆ k = 0 in general. The position of minimum, which isdenoted by ∆ k ∗ , is obtained by minimizing Eq. (A4) as∆ k ∗ = − Q c (cid:96) d (cid:96) u (cid:96) + Q d (cid:96) . (A5)The position of resonance shifts by this amount due tothe coupling Q .We denote by t (cid:0) x ± ( k ) , y ± ( k ) (cid:1) the two eigenvectorswith the corresponding eigenvalues (A3). The explicitforms of these eigenvectors are given by (cid:18) x ± ( k ) y ± ( k ) (cid:19) = 1 (cid:112) E (∆ k )[∆ E (∆ k ) ∓ a ] / (cid:18) ± ∆ E (∆ k ) − a ( − i) (cid:96) b (cid:19) . (A6)Thus the appropriate Floquet eigenstates are given by χ (cid:48) + ( k ) = x + ( k ) χ + ,(cid:96) ( k ) + y + ( k ) χ − , ( k ) , (A7) χ (cid:48)− ( k ) = x − ( k ) χ + ,(cid:96) ( k ) + y − ( k ) χ − , ( k ) (A8) instead of χ + ,(cid:96) ( k ) and χ − , ( k ) in the vicinity of reso-2nance.The new Floquet eigenstates χ (cid:48)± ( k ) carries harmoniccurrents, which are proportional to ∆ k . From Eqs. (19),(25), and (26), we obtain I resH ( k, N ) ≡ − ev N ( k ) (cid:88) α = ± | w α ( k ) | (cid:2) | x α ( k ) | − | y α ( k ) | (cid:3) (A9)= − ev N ( k ) u (cid:96) ∆ k ∆ E (∆ k ) (cid:2) | w − ( k ) | − | w + ( k ) | (cid:3) , (A10)where w ± ( k ) are now expansion coefficients in terms of χ (cid:48)± ( k ). This result reflects the fact that the Floqueteigenstates χ + ,(cid:96) ( k (cid:96) ) and χ − , ( k (cid:96) ) carry harmonic cur-rents with opposite sign.The HHC contribution near resonance I resH ( k, N ) be-comes as large as O ( Q ) if | u (cid:96) ∆ k | (cid:46) Q , although it van-ishes at an exceptional point ∆ k = 0, where | x ± ( k ) | = | y ± ( k ) | = 1 /
2. Thus, when summed over k , the HHCcontribution near resonance amounts to O ( Q ) owing tothe k -space volume factor of Q . We note that the k sumaround ∆ k does not vanish because ∆ E (∆ k ) is not aneven function of ∆ k and the denominator becomes min-imum at ∆ k = ∆ k ∗ (cid:54) = 0.The O ( Q ) contribution from resonance does not showa plateau for | F | (cid:46) N -dependence of I resH ( k, N ) derives from that of v N ( k ) andhence J N ( F ). As mentioned in Sec. IV, this does notshow a plateau for | F | (cid:46) O ( Q ). Mixing of χ (cid:48)± ( k ) with the other Floquet eigen-states with eigenvalues differring by M Ω ( M ∈ Z ) iscaused by Q at the first order, and this mixing leads to O ( Q ) contribution to the HHC spectrum. When summedover k , this contribution amounts to O ( Q ) due to the k -space volume factor of Q . One can show that the N -dependence of this contribution can be a plateau, but itswidth is about 2 z k rather than 4 z k obtained in Sec. V A.Thus the O ( Q ) contribution from resonances is not veryimportant to determine the width of the plateau in theHHC spectrum.Let us now numerically verify the Q -dependence of theHHC spectrum for small Q . For this purpose, we definethe HHC spectrum induced by the staggered potential∆ I H ( N ) ≡ I H ( N ) − I H ( N, Q = 0) , (A11)where I H ( N, Q = 0) denotes the HHC spectrum in thesingle-band limit discussed in Sec. IV. According to ourperturbation theory, this quantity is expanded as a poly-nomial in Q as∆ I H ( N ) = e ( a N Q + b N Q + · · · ) . (A12)We have shown that | a N | decreases faster than | b N | andthe O ( Q ) contribution becomes more important than -18 -14 -10 -6 -2 -4 -3 -2 (a) (b) Harmonic Order -15 -10 -5
1 3 5 7 9 11
FIG. 6. (a) Absolute value of ∆ I H ( N ) [Eq. (A11)] calculatedfor F = 0 . . Q in the log-logscale. Data are shown for the harmonics N = 1 (filled cir-cle), 3 (filled square), 5 (filled triangle), 7 (open circle), 9(open square), and 11 (open triangle). The solid lines arethe polynomial fit (A12) of the fourth order. (b) The de-termined fitting parameters normalized as a N /a (filled) and b N /b (open) plotted against the harmonic order N , where a = 6 . × − and b = − .
0. The error bars show thefitting errors. the O ( Q ) one for larger N . Figure 6(a) shows the numeri-cally calculated ∆ I H ( N ) for several Q ’s with F = 0 . .
0. The Q -dependence of ∆ I H ( N ) is fitted well bya polynomial (A12) of the fourth order, and this justifiesthe polynomial expansion. The first- and second-ordercoefficients a N and b N are shown in Fig. 6(b). This figureshows that the O ( Q ) contribution decreases with N fasterthan the O ( Q ) one, and this tendency is consistent withour analytical calculations. Thus these numerical datasupport our analysis with perturbation in Q includingthe resonance effects between the Floquet eigenstates. Appendix B: Second-Order Perturbation Theory
Here we extend the first-order perturbation theory inSec. V A to the second order. We derive the second or-der correction λ α, to the Floquet eigenstates, and showthat its O ( Q ) contribution to the HHC has a similar N -dependence to the one discussed in Sec. V A.The second-order correction λ α, ( k ) in Eq. (48) is asuperposition of χ α,M ( k ) with various M ’s with the same σ x -eigenvalue α since each V F mn [Eq. (50)] flips α . Thusthe correction is represented as λ α, ( k ) = (cid:88) M (cid:54) =0 C αM ( k ) χ α,M ( k ) , (B1)and the coefficients C αM ( k ) are given by the standard pro-cedure from the matrix elements (51) and the eigenener-gies as C αM ( k ) = ( α i) M M Ω (cid:88) M (cid:48) J M − M (cid:48) (2 z k ) J − M (cid:48) (2 z k ) M (cid:48) Ω + 2 αJ ( F ) cos k . (B2)One can easily check that C αM ( − k ) = ( − M C αM ( k ) = C − α − M ( k ).3We briefly interpret the M -dependence of C αM ( k ) inEq. (B2). We can safely ignore the vanishing of the de-nominator, which corresponds to the resonance discussedin Appendix A, because the contribution of the reso-nant region amounts to O ( Q ) due to the extra factor Q from the k -space volume. The sum in Eq. (B2) takesthe form of (cid:80) M (cid:48) g M (cid:48) J M − M (cid:48) (2 z k ) J − M (cid:48) (2 z k ) with a grad-ually changing g M (cid:48) . Now we recall that the n -dependenceof | J n (2 z k ) | is approximately constant for | n | (cid:46) z k anddecays rapidly for | n | (cid:38) z k . Therefore, the sum and,hence, C αM ( k ) depend on M rather slowly for | M | (cid:46) | z k | and rapidly decreases for | M | (cid:38) | z k | . We note thatthis behavior is not qualitatively modified by the overallfactor M − in Eq. (B2) since it varies slowly.Let us evaluate the second-order correction of the HHCspectrum from λ α, : I α, B H ( k, N ) ≡ (cid:88) n,l (cid:16) (cid:126)χ α, n − l + N ( k ) † I k,l (cid:126)λ α, n ( k )+ (cid:126)λ α, n − l + N ( k ) † I k,l (cid:126)χ α, n ( k ) (cid:17) . (B3)By invoking Eq. (B1) and performing some algebra, we obtain I α, B H ( k, N ) = − αe (cid:88) M (cid:54) =0 (cid:2) v N − M ( k ) C αM ( k ) + v ∗ N − M ( k ) C − αM ( k ) ∗ (cid:3) . (B4)Here v n ( k ) ∝ J n ( F ) has a significant weight only around n = 0 since | F | (cid:46)
1. Then the sum over M in Eq. (56)leaves αv ( k )[ C αN ( k )+ C − αN ( k ) ∗ ] and the N -dependence of I α, B H ( k, N ) is governed by that of C ± N ( k ). Since C ± N ( k )shows a plateau as shown above, I α, B H ( k, N ) also showsa plateau in its N -dependence. The width of the plateauis 4 | z k | .We remark that, apart from the resonance, the whole O ( Q ) contribution to the HHC spectrum is given bythe sum I α, ( k, N ) = I α, A H ( k, N ) + I α, B H ( k, N ), where I α, A H ( k, N ) derives from the first-order corrections to theFloquet eigenstates explained in Sec. V A and shows aplateau of width 4 | z k | . The analysis in this appendixshows that a similar plateau appears in the other part I α, B H ( k, N ) with the same width. [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] M. Hohenleutner, F. Langer, O. Schubert, M. Knorr,U. Huttner, S. W. Koch, M. Kira, and R. Huber, Nature , 572 (2015).[4] T. T. Luu, M. Garg, S. Y. Kruchinin, A. Moulet, M. T.Hassan, and E. Goulielmakis, Nature , 498 (2015).[5] G. Ndabashimiye, S. Ghimire, M. Wu, D. A. Browne,K. J. Schafer, M. B. Gaarde, and D. A. Reis, Nature , 520 (2016).[6] S. Ghimire, G. Ndabashimiye, A. D. DiChiara,E. Sistrunk, M. I. Stockman, P. Agostini, L. F. DiMauro,and D. A. Reis, Journal of Physics B: Atomic, Molecularand Optical Physics , 204030 (2014).[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] N. Yoshikawa, T. Tamaya, and K. Tanaka, Science (NewYork, N.Y.) , 736 (2017).[9] K. Kaneshima, Y. Shinohara, K. Takeuchi, N. Ishii,K. Imasaka, T. Kaji, S. Ashihara, K. L. Ishikawa, andJ. Itatani, Physical Review Letters , 243903 (2018).[10] T. Brabec and F. Krausz, Reviews of Modern Physics , 545 (2000).[11] D. Golde, T. Meier, and S. W. Koch, Journal of theOptical Society of America B , 2559 (2006).[12] D. Golde, T. Meier, and S. W. Koch, Physical ReviewB , 075330 (2008).[13] D. Golde, M. Kira, T. Meier, and S. W. Koch, PhysicaStatus Solidi B , 863 (2011).[14] G. Vampa, C. McDonald, G. Orlando, D. Klug, P. Corkum, and T. Brabec, Physical Review Letters ,073901 (2014).[15] T. Tamaya, A. Ishikawa, T. Ogawa, and K. Tanaka,Physical Review Letters , 016601 (2016).[16] T. Tamaya, A. Ishikawa, T. Ogawa, and K. Tanaka,Physical Review B , 241107 (2016).[17] L. Plaja and L. Roso-Franco, Physical Review B , 8334(1992).[18] M. Korbman, S. Yu Kruchinin, and V. S. Yakovlev, NewJournal of Physics , 013006 (2013).[19] M. Wu, S. Ghimire, D. A. Reis, K. J. Schafer, and M. B.Gaarde, Physical Review A , 043839 (2015).[20] T.-Y. Du and X.-B. Bian, Optics Express , 151 (2017).[21] T. Ikemachi, Y. Shinohara, T. Sato, J. Yumoto,M. Kuwata-Gonokami, and K. L. Ishikawa, Physical Re-view A , 043416 (2017).[22] G.-R. Jia, X.-H. Huang, and X.-B. Bian, Optics Express , 23654 (2017).[23] T. Ikemachi, Y. Shinohara, T. Sato, J. Yumoto,M. Kuwata-Gonokami, and K. L. Ishikawa,arXiv:1709.08153.[24] J. H. Shirley, Physical Review , B979 (1965).[25] N. Tzoar and J. I. Gersten, Physical Review B , 1132(1975).[26] F. H. M. Faisal and R. Genieser, Physics Letters A ,297 (1989).[27] N. B. Tal, N. Moiseyev, and A. Beswickf, Journal ofPhysics B: Atomic, Molecular and Optical Physics ,3017 (1993).[28] F. H. Faisal and J. Z. Kami´nski, Physical Review A ,R1769 (1996).[29] F. H. M. Faisal and J. Z. Kami´nski, Physical Review A , 748 (1997).[30] D. F. Martinez, L. E. Reichl, and A. L. A. Germ´an,Physical Review B , 174306 (2002). [31] J. Y. Yan, Physical Review B , 075204 (2008).[32] S. T. Park, Physical Review A , 013420 (2014).[33] T. N. Ikeda, Physical Review A , 063413 (2018).[34] A. K. Gupta, O. E. Alon, and N. Moiseyev, PhysicalReview B , 205101 (2003).[35] O. E. Alon, V. Averbukh, and N. Moiseyev, Advancesin Quantum Chemistry , 393 (2004).[36] H. Hsu and L. E. Reichl, Physical Review B , 1 (2006).[37] F. H. M. Faisal, J. Z. Kami´nski, and E. Saczuk, PhysicalReview A , 023412 (2005).[38] T. Higuchi, M. I. Stockman, and P. Hommelhoff, Phys-ical Review Letters , 213901 (2014).[39] D. Dimitrovski, T. G. Pedersen, and L. B. Madsen, Phys-ical Review A , 063420 (2017).[40] B. M. Breid, D. Witthaut, and H. J. Korsch, New Jour-nal of Physics , 110 (2006).[41] Y. Mizumoto and Y. Kayanuma, Physical Review A ,035601 (2012).[42] Y. Mizumoto and Y. Kayanuma, Physical Review A ,023611 (2013).[43] The energy gap minimum can be moved to k = 0 bya unitary transformation ˆ c j → ˆ c j exp(i π j ), which shiftsthe electron momentum k → k + π/ | Ψ (cid:105) = (cid:81) N el j =1 (cid:104) ˆ φ † k j · (cid:126)ψ k j (cid:105) | (cid:105) , where N el denotes the number ofelectrons and their momenta are { k , k , . . . , k N el } .[46] We implicitly assume A ( t ) = A sin Ω t for −∞ < t < ∞ in using the Floquet theory. Once we impose the initialcondition (cid:126)ψ k (0) = (cid:126)ψ k , the solution of Eq. (11) is appro-priately obtained in t ≥ t ≥ I ( t <
0) = 0.[48] To prove this, we note the following sum rule: I (0) = (cid:80) N I H ( N ) + (cid:82) ∞−∞ d ω π I C ( ω ). Here, (cid:82) ∞−∞ d ω π I C ( ω ) and,hence (cid:82) ∞−∞ d ω π | I C ( ω ) | are finite since both I (0) and (cid:80) N I H ( N ) are finite. Thus it follows from the Parseval-Plancherel identity (cid:82) ∞−∞ d t | I C ( t ) | = (cid:82) ∞−∞ d ω π | I C ( ω ) | that lim t →∞ I C ( t ) = 0.[49] K. A. Pronin, A. D. Bandrauk, and A. A. Ovchinnikov,Physical Review B , 3473 (1994).[50] Near k = 0, rapid decays actually occur for | N | (cid:46) | z k | .This is because M Ω plays a minor role in Eq. (53) com-pared with 2 αJ ( F ) cos k . In fact, if we neglect M Ω, wehave I α, ( k, N ) = − α [2 J ( F ) cos k ] − v N ( k ) ∝ J N ( F ),which rapidly decays with N . Thus the plateau of theHHC spectrum is contributed from larger k values.[51] K. Shimomura, K. Uchida, K. Nagai, and K. Tanaka (un-published).[52] R. E. Silva, I. V. Blinov, A. N. Rubtsov, O. Smirnova,and M. Ivanov, Nature Photonics , 266 (2018),arXiv:1704.08471.[53] Y. Murakami, M. Eckstein, and P. Werner,arXiv:1712.06460 (2017).[54] Y. Murakami and P. Werner, arXiv:1804.08257 , 1 (2018).[55] T. Nag, R.-J. Slager, T. Higuchi, and T. Oka,arXiv:1802.02161 (2018).[56] K. E. Hamilton, A. A. Kovalev, A. De, and L. P.Pryadko, Journal of Applied Physics (2015),10.1063/1.4921929.[57] Y. Zhou and M. W. Wu, (2011), 10.1103/Phys-RevB.83.245436, arXiv:1103.4704.[58] T. Morimoto and N. Nagaosa, Science Advances ,e1501524 (2016).[59] T. Morimoto and N. Nagaosa, Physical Review B ,035117 (2016).[60] K. W. Kim, T. Morimoto, and N. Nagaosa, PhysicalReview B , 035134 (2017).[61] H. Yamane and S. Tanaka, Symmetry10