High-order harmonic generation in graphene: nonlinear coupling of intra and interband transitions

We investigate high-order harmonic generation (HHG) in graphene with a quantum master equation approach. The simulations reproduce the observed enhancement in HHG in graphene under elliptically polarized light [N. Yoshikawa et al, Science 356, 736 (2017)]. On the basis of a microscopic decomposition of the emitted high-order harmonics, we find that the enhancement in HHG originates from an intricate nonlinear coupling between the intraband and interband transitions that are respectively induced by perpendicular electric field components of the elliptically polarized light. Furthermore, we reveal that contributions from different excitation channels destructively interfere with each other. This finding suggests a path to potentially enhance the HHG by blocking a part of the channels and canceling the destructive interference through band-gap or chemical potential manipulation.

High-order harmonic generation (HHG) is an extreme photon-upconversion process based on highly nonlinear light-matter interactions. HHG was originally observed in atomic gas systems more than thirty years ago [1,2]. Several years later, the microscopic mechanism underlying HHG in noble gases was beautifully explained by a simple semi-classical model, the so-called three-step model [3][4][5]. On the practical side, high coherence in these upconversion processes allows us to generate extremely short light pulses, presenting a novel avenue to time-domain investigations of ultrafast electron dynamics in matter [6][7][8][9][10][11][12][13][14].
Since the first observation of HHG in solids by Gimire et al. [15], HHG in solids has been attracting much interest as it may have various applications ranging from the development of novel light sources [16] to probing of microscopic information of matter [17][18][19]. So far, experimental studies on HHG in solids have explored various materials [20][21][22][23][24], and the theoretical aspects of HHG in solids have been intensively investigated with various approaches [25][26][27][28][29][30].
Recently, attosecond transient absorption spectroscopy and microscopic simulations clarified that nonlinear coupling of intraband and interband transitions play significant roles in ultrafast modification of optical properties and in nonlinear photocarrier-injection processes [11,31,32]. Hence the nonlinear coupling of the two kinds of transitions may be a key to accessing the microscopic physics behind a light-induced phenomenon and may offer a novel opportunity to control it. While intraband and interband transitions have been discussed in the context of HHG in solids [25,33], still detailed roles of nonlinear coupling of these transitions have not yet been investigated. Yoshikawa et al. recently reported that the HHG in graphene can be enhanced by elliptically-polarized light [23]. This observation is distinct from the HHG in noble gases, where HHG is significantly suppressed with an increase in the ellipticity of light [34][35][36]. Therefore, HHG in graphene under elliptically-polarized light would offer an opportunity to look into the microscopic mechanism underlying HHG in solids. However, the microscopic mechanism of the HHG enhancement under ellipticallypolarized light has been still unclear.
In this work, we investigate the enhancement of HHG in graphene with elliptically-polarized light, by employing a quantum master equation with a simple two-band model. The simple modeling of electron dynamics in graphene fairly captures the experimentally-observed enhancement of HHG and provides a microscopic insight into the mechanism. The model indicates a significant role of the nonlinear coupling between light-induced intraband and interband transitions in the enhancement of HHG, demonstrating a destructive interference among multiple HHG channels.
To describe the electron dynamics, we employ a quantum master equation with a two-band approximation for the Dirac cone of graphene [37,38]. In the model, the time propagation of the one-body reduced density matrix at each Bloch wavevector k is described by where H k+A(t) is the Hamiltonian, andD[ρ k (t)] is a relaxation operator. In this work, we employ the following 2-by-2 Hamiltonian matrix: arXiv:2010.06275v1 [cond-mat.mtrl-sci] 13 Oct 2020 where σ j are Pauli matrices, k j are the j-component of the Bloch wavevector k, and A j (t) is the j-component of the vector potential A(t), which corresponds to the applied electric fields, E(t) = −Ȧ(t). The band-gap ∆ is set to zero for graphene unless stated otherwise. Here, τ z determines the chirality of the system (either +1 or −1). We evaluate observables as the average of two calculations with opposite chiralities. We set the Fermi velocity v F to 1.12 × 10 6 m/s in accordance with an abinitio simulation [39]. The relaxation operatorD[ρ k (t)] is constructed by making the relaxation-time approximation with a longitudinal relaxation time of T 1 = 100 fs and transverse relaxation time of T 2 = 20 fs [37,38,40]. Note that the relaxation operator also depends on the chemical potential µ; µ is set to zero (charge neutrality point) unless stated otherwise.
To describe the applied electric fields, we employ the following form of the vector potentials in the domain −T f ull /2 < t < T f ull /2; the potential is zero outside this domain. In accordance with the experimental conditions in Ref. [23], we set the mean photon energy ω 0 to 260 meV. The full pulse duration T f ull is set to 100 fs. We will investigate the electron dynamics by changing the peak field strength of the applied laser fields, E 0,x and E 0,y . Employing a time-dependent density matrix, ρ k (t), we can evaluate the induced electric current as whereĴ is the current operator defined bŷ Note that the current defined in Eq. (4) depends on E 0,x and E 0,y via A(t) in Eq. (3), and for clarity we will indicate this dependence in the next equation by using the notation J (t, E 0,x , E 0,y ). In experiments, HHG occurs not only at the center of the beam-spot but also on the whole focal area. To make our model more realistic, we employ the following intensity-averaging procedure to approximate the results for the case of a Gaussian beam profile [28,40]: The power spectrum of the high-order harmonics polarized along the j-direction can be evaluated with the current as where J ave j (t) is the j-component of the current vector J ave (t). Furthermore, the intensity of the nth-order harmonics can be evaluated by integrating the power spectrum within a finite range as First, we evaluate the ellipticity dependence of the HHG by fixing the peak field strength E 2 0,x + E 2 0,y to 6.5 MV/cm inside the material. The major axis of the elliptically polarized light is set to the x-axis while the minor axis is set to the y-axis. Figure 1 (a) shows the signal intensity of the 7th-order harmonics as a function of laser ellipticity, E 0,y /E 0,x . The intensity I 7th j is separately computed for the different polarization directions (j = x or j = y) of the emitted harmonics. As seen from the figure, when the applied laser field is linearly polarized in the x-direction, the emitted high harmonics are also linearly polarized in the x-direction. Once the applied fields become elliptical, the emitted harmonics also become elliptical, having both x and y-components. Interestingly, the y-component I 7th y rapidly increases with the increase in driver ellipticity and becomes much larger than the x-component I 7th x , demonstrating enhancement in HHG in graphene by elliptically-polarized light. Once the driver ellipticity further increases and approaches one (circularly-polarized light), the emitted harmonics is significantly suppressed due to the circular symmetry of the Dirac cone. These observations are consistent with the experimental results [23]. Thus, it has been demonstrated that the simple Dirac band model with the relaxation-time approximation contains sufficient ingredients to describe the HHG in graphene and can be used to investigate its microscopic origin.
To obtain further insight into the phenomena, we evaluated the harmonic intensity by changing the band-gap ∆. Figure 1 (b) shows the 7th-order harmonic intensity as a function of band-gap ∆. Here, we used the same field strength as in Fig. 1 (a). The ellipticity is set to 0.17, by the harmonic intensity is maximized at ∆ = 0. Surprisingly, the harmonic intensity can be significantly enhanced by increasing the band-gap. Furthermore, the x-component of the harmonic intensity I 7th x shows a peak around a band-gap of 0.8 eV, which is close to the energy of three photons (0.78 eV), while the y-component I 7th y shows a peak around a band-gap of 0.5 eV, which is close to the energy of two photons (0.52 eV). The enhancement and formation of peaks that occur as the band-gap increases indicate that multi-photon processes play a significant role in HHG in graphene, while Zener tunneling is expected to have only a minor contribution in the present regime.
The enhancement in HHG with the increase in bandgap ∆ may be regarded as a counter-intuitive consequence because an increase in the gap tends to block a The 7th-order harmonic intensity from graphene under elliptically polarized light: (a) harmonic intensity as a function of ellipticity, E0,y/E0,x. The experimental data [23] are also shown: (b) harmonic intensity as a function of the band-gap ∆. The harmonic intensity is decomposed into the two polarization axes: the major I 7th x (red-solid) and minor I 7th y (blue-dotted) axes of the driving elliptically polarized light.
part of the transitions. In fact, the high-order harmonics vanish once the band-gap becomes significantly large, as shown in Fig. 1 (b). To understand the enhancement in HHG with the increase in the gap, we propose a microscopic mechanism based on destructive interference between multiple channels: high-order harmonics are generated as a superposition of multiple signals from various microscopic paths due to nonlinear coupling of intraband and interband transitions. We further suppose that the multiple signals may destructively interfere with each other, and the total signal may be weakened. When such destructive interference plays a significant role, HHG may be enhanced by increasing the gap, because in so doing contributions can be partly suppressed, and the destructive interference can be canceled.
To examine our hypothesis, let us investigate the interference of different HHG contributions from the viewpoint of intraband and interband transitions. Here, we partly turn off the transitions based on the instantaneous eigenbasis representation, where the intraband transitions appear in the diagonal elements of the Hamiltonian and the interband transitions appear in the off-diagonal elements [40,41]. Elliptically-polarized light consists of two polarization components, and each of them induces intraband and interband transitions. Hence, we can con-sider the four kinds of light-induced transitions. For later convenience, we label them as follows: intraband transitions induced by the x-component of the electric fields (τ a ) and by the y-component (τ b ); likewise, interband transitions induced by the x-component of the electric fields (τ c ) and by the y-component (τ d ). Figures 2 (a-d) show the strength distributions of each transition in kspace. Here, the strength of the intraband transitions is evaluated as the gradient of the single-particle energy, To elucidate the roles of the intraband and interband transitions in HHG, we can compute the electron dynamics by turning off part of them. In Fig. 2 (e), the 7th-order harmonic intensity for the y-direction, I 7th y , with both intraband and interband transitions is shown as the blacksolid line. The results for only intraband transitions are shown as the red dashed line, while those for only interband transitions are shown as the blue dotted line. Here, we have used the same conditions as in Fig. 1 (b) except the band-gap ∆, which is set to a small value of 0.035 eV to avoid numerical singularity in the intrabandinterband transition analysis at the Dirac point. As the figure clearly shows, neither pure intraband nor interband transitions can induce the HHG. Therefore, HHG originates from a nonlinear coupling of interband and interband transitions.
To study the roles of the coupling among the intraband and interband transitions, we can decompose the current J ave (t) into the coupling components of the transitions as follows: where τ , σ, and δ denote the labels of the transitions (τ a , τ b , τ c , and τ d ). In Eq. (9), there are 15 terms with four kinds of current: J ave τ (t) is current induced solely by the transition τ . J ave τ,σ (t) is current induced by the coupling of two transitions, τ and σ. Likewise, J ave τ,σ,δ (t) is the current induced by the coupling among three transitions, τ , σ, and δ. Finally, J ave τa,τ b ,τc,τ d (t) is the current induced by the coupling of all four transitions. For more details, see Supplemental Material [40].
We evaluated the high-order harmonic intensity with the decomposed currents in Eq. (9) instead of the total current J ave (t). Figure 2 (f) shows the 7th-order harmonic intensity I 7th y as a function of ellipticity for various decomposed current. Here, only the five major contributions, J ave τa,τ d (t), J ave τa,τ b ,τc (t), J ave τa,τ b ,τ d (t), J ave τa,τc,τ d (t), J ave τa,τ b ,τc (t) and J ave τa,τ b ,τc,τ d (t), are shown, while all the other contributions are rather minor [40]. In fact, the reconstructed signal from the five major contributions (grey dotted line) shows fair agreement with the full signal (black solid line). Remarkably, all the decomposed re-sults in Fig. 2 (b) have a larger harmonic intensity than the total signal. Hence the results clearly demonstrate destructive interference among the various contributions. Therefore, our hypothesis, i.e., destructive interference of HHG, is clearly supported by the theoretical results.
As seen from Fig. 2 (f), the coupling of the intraband transition induced by the x-component of the electric fields (τ a ) and the intraband transitions induced by the y-component of the fields (τ d ) shows the largest contribution to the harmonic intensity (red-dashed line). Therefore, the cross-coupling between the intraband and interband transitions induced respectively by the perpendicular components of the electric fields plays a significant role in the enhancement of HHG in graphene under elliptically-polarized light. In fact, all five major contributions include cross-coupling of the intraband and interband transitions with the perpendicular field components. This observation can be straightforwardly understood in terms of the transition strength distribution in Figs Having established the destructive interference mechanism of HHG in solids, we propose a way to enhance HHG by canceling the destructive interference with the chemical potential shift. Since the chemical potential shift suppresses a part of the transitions in graphene by Pauli blocking (see the inset of Fig. 3), the destructive interference of multiple channels may be canceled by the tuning of the chemical potential, and this should result in an enhancement of HHG. To demonstrate the impact of the chemical potential shift, we evaluate the 7th-order harmonic intensity from graphene by changing the chemical potential µ. Figure 3 shows the 7th-order harmonic intensity as a function of the chemical potential µ. Here, we used the same field strength as Fig. 1 (b) and set the ellipticity to 0.17. As shown in Fig. 3, the highorder harmonic intensity can be enhanced by tuning the chemical potential. Hence the proposed method of enhancement in the HHG based on the destructive interference mechanism has clearly been demonstrated. Furthermore, the harmonic intensity shows the peak around |µ| ∼ 1 eV ∼ 7 ω/2. Note that, once the absolute value of the chemical potential reaches half of the n-th-photon energy, |µ| ∼ n ω 0 /2, the n-photon resonant processes are suppressed by the Pauli blocking. Thus, the enhancement in HHG and the peak feature in Fig. 3 further indicate the significance of multi-photon processes in the present regime. In conclusion, we investigated and identified the microscopic mechanism underlying HHG in graphene exposed to elliptically polarized light by employing the quantum master equation with the simple Dirac cone [37,40]. We found that the nonlinear coupling between the intraband and interband transitions is the microscopic origin of the HHG in graphene. In particular, the cross-coupling between the intraband and interband transitions induced by the perpendicular components of the electric fields causes the enhancement in HHG under elliptically polarized light, reflecting the unique transition strength profiles of graphene in k-space, as shown in Figs. 2 (a-d).
Our findings of the interference effects on HHG will lead to a general understanding of HHG in solids, and also provide novel techniques to increase up the HHG efficiency. Indeed, we have demonstrated that the highorder harmonics can be enhanced by tuning the chemical potential µ or the band-gap ∆, by blocking a part of the transitions. Furthermore, the interference mechanism may open a way to control the high-order harmonic generation with a large degree of freedom by tuning phases of the HHG channel contributions and by flipping the destructive interference to constructive interference through optimization of the applied laser fields. where T 1 and T 2 are the time-constants for the longitudinal and transverse relaxations, respectively. In this work, as a target state of the relaxation, we employ the Fermi-Dirac distribution where µ is the chemical potential, T e is the electron temperature. In this work, we fix the electron temperature T e to 80 K.

II. FOCAL SPOT AVERAGE
In realistic experimental configurations, the high-order harmonic generation occurs not only at the center of a beam-spot but also on the whole focal area. Thus, the generated high-order harmonics from a wide region of the focal area can be detected. To take into account the macroscopic focal-spot average effect of HHG, we employ the intensity average procedure according to Ref. [6].
Here, we assume the field strength of laser electric fields has the following Gaussian profile on the sample surface where (x, y) are the coordinates on the surface, E 0 is the peak field strength, and σ is the beam waist. We further assume that the beam waist is sufficiently large, and the induced current depends only on the local field strength as J [E(x, y)](t). Based on these assumptions, the average current within the beam spot can be evaluated as Hence the total current on the sample can be evaluated as the intensity average of the induced current. In this work, we repeatedly perform the simulation by changing the laser field strength E 0 = E 2 0,x + E 2 0,y with fixed ratio of E 0,x and E 0,y . Then, we compute the averaged current on the sample J ave (t) with Eq. (9).
band. Hence the transitions induced by A tra (t) is nothing but the intraband transitions. On the other hand, if A tra (t) = 0 and A ter (t) = A(t), the system shows only a transition among different bands without the motion in the k-space. Hence the transitions induced by A ter (t) is nothing but the interband transitions. Note that these definitions of the intraband and interband transitions with the reduced density matrix are equivalent to those defined with the Houston expansion in the wavefunction theory [7].
Based on the electron dynamics simulations with the effective HamiltonianH tra−ter k+A(t) , one can elucidate impacts of intraband and interband transitions. For example, by setting A tra (t) to A(t) and A ter (t) to zero, one can study the electron dynamics induced solely by the intraband transitions. Likewise, by setting A tra (t) to zero and A ter (t) to A(t), one can study the electron dynamics induced solely by the interband transitions. Applying this approach to the HHG in graphene, we evaluated the high-order harmonic intensity induced solely by intraband or interband transitions. The results are shown in Fig. 2 (e) in the main text.
To investigate detailed roles of nonlinear coupling among intraband and interband transitions, we introduce a decomposition of the current J (t, E 0,x , E 0,y ) into each nonlinear coupling component. In the above analysis, we introduced the two gauge fields, A tra (t) and A ter (t), to distinguish the intraband and interband transitions. Furthermore, each gauge field vector consists of two directional components. Hence, the induced current is a function of the field strength of each component of each gauge field as J tra−ter (t, E 0,x,tra , E 0,x,ter , E 0,y,tra , E 0,y,ter ). In the main text, we introduced the following four labels for each kind of transitions. τ a : the intraband transitions induced by the xcomponent of fields. τ b : the interband transitions induced by the x-component of fields. τ c : the intraband transitions induced by the y-component of fields. τ d : the interband transitions induced by the y-component of fields. With this notation, we first define four kinds of current as follows: J τc (t) = J tra−ter (t, 0, 0, E 0,y,tra , 0), J τ d (t) = J tra−ter (t, 0, 0, 0, E 0,y,ter ).
We evaluated the harmonic intensity with each decomposed current. Figures S1-S3 show the 7th-order harmonic intensity I 7th y as a function of ellipticity for various decomposed current. Here, we employed the same conditions as those of Fig. 2 (f) of the main text. In Fig. S1, the results of the current induced sorely by a single transition, J τ , are shown. In Fig. S2, the results of the current induced by the nonlinear coupling between two of four kinds of transitions, J τ σ , are shown. In Fig. S3, the results of the current induced by the nonlinear coupling among three of four transitions, J τ σδ , and among all of four transitions, J τa,τ b ,τc,τ d (t), are shown. In these figures, the result of the full transitions is also shown as the black-solid line.