Enhancement of high-order harmonic generation in graphene by mid-infrared and terahertz fields

We theoretically investigate high-order harmonic generation (HHG) in graphene under mid-infrared (MIR) and terahertz (THz) fields based on a quantum master equation. Numerical simulations show that MIR-induced HHG in graphene can be enhanced by a factor of 10 for fifth harmonic and a factor of 25 for seventh harmonic under a THz field with a peak strength of 0.5 MV/cm by optimizing the relative angle between the MIR and THz fields. To identify the origin of this enhancement, we compare the fully dynamical calculations with a simple thermodynamic model and a nonequilibrium population model. The analysis shows that the enhancement of the high-order harmonics mainly results from a coherent coupling between MIR- and THz-induced transitions that goes beyond a simple THz-induced population contribution.

Among various materials, the HHG from graphene has been intensively investigated theoretically.It has been suggested that HHG can be efficiently induced in graphene because of the unique electronic structure of this material, that is, the Dirac cone [24][25][26][27].Recently, HHG from graphene has been experimentally observed in the mid-infrared (MIR) [28,29] and terahertz (THz) [30,31] regimes, exhibiting unique ellipticity dependence and high efficiency.We previously investigated HHG from graphene in both the MIR and THz regimes based on a quantum master equation [32,33].In the MIR regime, coupling between field-induced intraband and interband transitions opens important channels for HHG, enhancing HHG with finite ellipticity [32].A realtime electron dynamics simulation in the THz regime has shown that it is essential to consider the nonequilibrium steady-state resulting from the balance between * ssato@ccs.tsukuba.ac.jp field-driving and relaxation to go beyond the equilibrium thermodynamic picture of HHG from graphene [33].
It is important to improve the efficiency of solid-state HHG to develop novel HHG-based light sources and spectroscopies.In recent studies, it has been suggested that HHG from graphene can be enhanced using twocolor laser fields based on various mechanisms [34][35][36].In Ref. [34], the two-color HHG is suggested with the combination of the electron-hole pair creation by highfrequency pump light and the acceleration of the created pairs by low-frequency light.Mrudul et al investigated the HHG from graphene with bicircular fields, controlling the valley polarization [35].Furthermore, Avetissian et al investigated the HHG from graphene with a linearly polarized light and its second harmonics, showing that when the two-color fields are perpendicularly polarized, the stronger harmonics can be emitted compared with the parallel polarization [36].
In this study, we explore the possibility of using a THz field to enhance MIR-induced HHG in graphene based on the knowledge gained from previous studies.First we use a quantum master equation to compute the electron dynamics under MIR and THz fields and evaluate the emitted harmonic spectra.We compare the results of the fully dynamical calculations with a thermodynamic model considering the equilibrium Fermi-Dirac distribution and with a nonequilibrium population model considering a population distribution in a nonequilibrium steady-state.We find that coherent driving by the THz field plays an essential role in enhancing MIR-induced HHG beyond the induced population effect.
The paper is organized as follows.In Sec.II, we briefly describe the theoretical method for describing electron dynamics in graphene induced by MIR and THz fields.In Sec.III, we investigate the impact of a THz field on MIR-induced HHG from graphene and explore the microscopic mechanism of the enhancement by employing the thermodynamic model and the nonequilibrium population model.Finally, our findings are summarized in Sec.IV.

A. Electron dynamics calculation based on a quantum master equation
Here, we briefly describe theoretical methods for computing light-induced electron dynamics in graphene and how the calculated dynamics can be used to analyze HHG.These methods have been described in detail in previous works [33,37,38].
In this study, we use the following quantum master equation to describe electron dynamics in graphene: where ρ k (t) is the one-body reduced density matrix for the Bloch wavevector k, and H k+eA(t)/ is the one-body Hamiltonian.Here, light-matter coupling is described by an additional term in the Hamiltonian, a spatially uniform vector potential A(t), which is related to the applied electric field as A(t) = − t −∞ dt E(t ), in the long-wavelength approximation.The relaxation operator is denoted as D [ρ k (t))].Based on the previous studies [33,37], we employ the relaxation time approximation [39] for the relaxation operator in the Houston basis expression [40,41].Those Houston states |u H bk (t) are simply the instantaneous eigenstates of the timedependent Hamiltonian and satisfy the following instantaneous eigenvalue problem: where b denotes the band index, and b,k+eA(t)/ are the instantaneous eigenvalues of the Hamiltonian.The reduced density matrix can be expanded in the Houston states as where ρ bb ,k (t) are the expansion coefficients.The relaxation operator D [ρ k (t)] can also be expanded in the Houston states as From Refs.[33,37], the longitudinal relaxation time T 1 is set to 100 fs, and the transverse relaxation time T 2 is set to 20 fs.Here, f FD ( , T e , µ) = 1 e ( −µ)/k B Te +1 is the Fermi-Dirac distribution, in which we set the electron temperature T e to 300 K and the chemical potential µ to 0, unless stated otherwise.
We describe the electronic structure of graphene employing a tight-binding Hamiltonian with nearestneighbor hopping [42] as follows: where t 0 is the nearest-neighbor hopping, and f (k) is given by f (k) = e ik•δ1 + e ik•δ2 + e ik•δ3 with the nearestneighbor vector δ j .According to Ref. [[42]], the nearestneighbor hopping t 0 is set to 2.8 eV, and the lattice constant a is set to 1.42 Å.By using the time-dependent density matrix ρ k (t) evolved with Eq. (1), the total energy of the electronic system can be evaluated as where Ω BZ is the volume of the first Brillouin zone.
Similarly, the electric current density is given by where the current operator Ĵ k (t) is defined as The current density J (t) induced by an intense electric field E(t) is analyzed to investigate HHG.For example, the power spectrum of the emitted harmonics can be evaluated by applying the Fourier transform to the current density J (t) as follows: III. RESULTS

A. MIR-induced HHG in graphene under THz fields and the quasistatic approximation
We first analyze the HHG induced by a MIR laser pulse in the presence of THz fields.For practical calculations, we employ the following form for the MIR pulse: in the domain −T MIR /2 < t < T MIR /2, and zero outside this domain.Here, E MIR is the peak strength of the MIR field, ω MIR is the mean frequency, e MIR is a unit vector along the polarization direction of light, and T MIR is the pulse duration.In this study, the pulse duration T MIR is set to 40 ps, and the mean frequency ω MIR is set to 0.35424 eV/ .We compute the electron dynamics by changing the other parameters.Similarly, we employ the following form for the THz pulse, in the domain −T THz /2 < t < T THz /2, and zero outside this domain.
Here, E THz is the peak strength of the THz field, ω THz is the mean frequency, e THz is a unit vector along the polarization direction, and T THz is the pulse duration.In this study, the pulse duration T THz is set to 42 ps and the mean frequency ω THz is set to 1.2407 meV/ .The time profile of the applied THz electric field is shown in the inset of Fig. 1 (a).
To gain insight into THz-assisted MIR-induced HHG in graphene, we perform the electron dynamics calculation in the presence of both THz and MIR fields, E THz (t) + E MIR (t).Here, we set E MIR to 6.5 MV/cm and E THz to 0.5 MV/cm.We note that intense THz pulses with amplitudes exceeding 1 MV/cm are available experimentally [43].The polarization direction of the THz field e THz is set to the Γ-M direction (the xdirection in our setup), whereas the polarization direction of the MIR field e MIR is treated as a tunable parameter.Figures 1 (a We analyze MIR-induced HHG by extracting the current induced by the MIR field in the presence of the THz field.We define two kinds of currents for this purpose.The first current is induced by both the THz and MIR fields and is denoted as J THz+MIR (t).The second current is induced solely by the THz field and is denoted as J THz (t).We define the current induced by the MIR field in the presence of the THz field as J eff (t) = J THz+MIR (t) − J THz (t).The Fourier transform is applied to the extracted current J eff (t), and the power spectrum of the emitted harmonics is computed using Eq. ( 9).The dashed line in Fig. 1 (e) shows the power spectrum computed using the current J (t) presented in Fig. 1 (a), where the polarization directions of the THz and MIR fields are parallel, and the dashed line in Fig. 1 (f) shows the power spectrum computed using the current J (t) presented in Fig. 1 (b), where the polarization directions of the two fields are perpendicular.Figure 1 (e) shows that second and higher even-order harmonics are generated in addition to odd-order harmonics because the THz field breaks the inversion symmetry of the system locally in time.This second-harmonic generation is known as electric-field-induced second-harmonic generation (EFISH) or THz-induced second-harmonic generation (TFISH) [44][45][46][47].Even-order harmonics are similarly generated in the perpendicular configuration (e MIR ⊥ e THz ), as can be seen from Fig. 1 (f).
Explicit use of the THz pulse in the electron dynamics calculation increases the propagation time (4000 fs in the present case), as can be seen Figs. 1 (a) and (b).Hence, the electron dynamics calculation with the explicit inclusion of THz pulses has a large computational cost.To lower the computational cost of modeling MIR-induced HHG in graphene in the presence of a THz field, we replace the THz pulse with a static field to calculate the electron dynamics, corresponding to a quasistatic approximation [33].For practical analysis, we perform two simulations.In the first simulation, the electron dynamics are calculated under a static field E dc (t) = e dc E dc Θ(t) that is suddenly switched on at t = 0. Here, e dc is the unit vector along the polarization direction of the static field, and E dc is the field strength.Immediately after the static field is switched on, the driven electron dynamics induce a current.The driven system reaches a steady state after a sufficiently long time propagation time, and the current becomes constant in time.We denote the current induced solely by E dc (t) as J dc (t).In the second simulation, the electron dynamics are calculated under both MIR and static fields, E dc (t) + E MIR (t − τ MIR ), where the pulse center of the MIR field is shifted by τ MIR .We denote the current induced by E dc (t) + E MIR (t − τ MIR ) as J dc+MIR (t).The shift τ MIR can be made sufficiently large time to investigate the MIR-induced electron dynamics for a full nonequilibrium steady state realized by the static field E dc (t).The MIR-induced current can be extracted as J eff (t) = J dc+MIR (t) − J dc (t) to analyze MIR-induced HHG in the presence of the static field.
Figures 1 (c) and (d) show the current J dc+MIR (t) induced by both static and MIR fields.Here, the static field is polarized along the Γ-M direction (the x-direction in our setup), and the field strength E dc is the same as the peak strength of the THz field, E dc = E THz = 0.5 MV/cm.The MIR field used in Fig. 1 (c) is the same as that used in Fig. 1 (a) and has a polarization direction parallel to that of the static field.By contrast, the MIR field used in Fig. 1 (d) is the same as that used in Fig. 1 (b) and has a polarization direction perpendicular to that of the static field.To apply the MIR field to the nonequilibrium steady-state under the static field, we set the time delay τ MIR of the MIR field to 1 ps, which is sufficiently longer than the relaxation time scales of the (1) (1) quantum master equation, T 1 and T 2 .
To analyze HHG in the presence of the static field E dc (t), we extract the current J eff (t) induced by the MIR field in the presence of the static field by subtracting J dc (t) from J dc+MIR (t): J eff (t) = J dc+MIR (t) − J dc (t).The solid lines in Figs 1 (e) and (f) correspond to the HHG spectra computed using the current shown in Figs 1 (c) and (d), respectively.The results of the quasistatic approximation are identical to those computed by explicitly including the THz pulse.Therefore, the quasistatic approximation is very good for analyzing HHG under MIR and THz fields.Hereafter, we employ the static field within the quasistatic approximation instead of explicitly including the THz pulse.The agreement between results obtained using the quasistatic approximation and the explicit inclusion of the THz pulse indicates that the nonequilibrium steady state under the static field plays an important role in MIR-induced HHG in graphene in the presence of a THz field.

B. Orientational dependence of HHG
Here, we investigate HHG in graphene within the quasistatic approximation by changing the relative angle between the static and MIR fields.For practical analysis, the direction of the static field e dc is fixed to the Γ-M axis (the x−axis in our setup), and the peak field strength of the MIR field E MIR is fixed at 6.5 MV/cm.The emitted harmonics are investigated by manipulating the polarization direction of the MIR field, e MIR , and the strength of the static field, E dc .
To analyze the HHG efficiency, we compute the signal intensity of the emitted harmonics at each order by integrating the power spectrum within a finite energy range as follows: Here, I nth total is the integrated intensity of the emitted nth harmonic.
Figures 2 (a-d) show the computed angular dependence of the emitted harmonic yield I nth for different harmonic orders.The angle θ denotes the relative angle between the MIR and static fields.In Fig. 2 (a), there is no second harmonic because graphene has intrinsic inversion symmetry in the absence of a static field.By contrast, the second harmonics are generated under the application of a static field because of the breakdown of the inversion symmetry.For a static field strength of 0.5 MV/cm, the emitted second-harmonic intensity is maximized at a relative angle of approximately 45 • .
In Fig. 2 (b), the third-harmonic yield is almost isotropic (black line) in the absence of a static field, reflecting the rotational symmetry of the Dirac cone (see also Appendix A).By contrast, the third-harmonic intensity exhibits a strong angular dependence under the application of a strong static field (E dc = 1.0 MV/cm): the third-harmonic emission is considerably enhanced when the static and MIR fields are perpendicular to each other and suppressed when these two fields are parallel.The enhancement of the third harmonic for the perpendicular configuration can be understood in terms of the coupling between the intraband transition induced by the static field and the interband transition induced by the MIR field, as was suggested in a previous study [32].
By comparison, the higher-order harmonics exhibit a more complex angular dependence under a static field, as shown in Fig. 2 (c and d).In particular, the fifth-order harmonic emission can be considerably enhanced in the presence of static or THz fields (Fig. 2 (d)).For example, the intensity of the fifth-order harmonic is enhanced more than ten times by applying a static field with a strength of 0.5 MV/cm with respect to the result solely induced by the MIR field (see the green line in Fig. 2 (d)).This enhancement ratio is larger than that of the thirdorder harmonic.Hence, a larger field-induced enhancement is expected for higher-order harmonics.In fact, the seventh-order harmonic is enhanced 25 times when the static field strength is 0.5 MV/cm (see Appendix B).
To further elucidate the angular dependence of HHG in graphene, we decompose the harmonic intensity I HHG (ω) into parallel and perpendicular components with respect to the polarization of the driving MIR field.The parallel component of the HHG intensity is defined as where e MIR is the unit vector along the polarization direction of the MIR field.Likewise, the perpendicular component is defined as where ēMIR is a unit vector perpendicular to e MIR , i.e., ēMIR • e MIR = 0.The total intensity I HHG in Eq. ( 9) is reproduced by the sum of I para HHG (ω) and I perp HHG (ω) as I HHG (ω) = I para HHG (ω) + I perp HHG (ω).Eq. ( 13) and Eq. ( 14) are directionally decomposed to analyze the parallel and perpendicular components of the emitted harmonics with respect to the polarization direction of the MIR field.Figures 2 (e-h) and (i-l) show the angular dependence of the parallel and perpendicular components the harmonic intensity for orders, respectively.
In Figs 2 (a), (e), and (i), the parallel component of the second harmonic under the static field reaches a maximum at approximately 45 • and dominates the total second-harmonic intensity at this angle.By contrast, the maximum perpendicular component is always obtained when the MIR and static fields are perpendicular to each other.In Figs. 2 (b), (f), and (j), the third harmonic is dominated by the parallel component for any angle and static field strength over the investigated range.For both second-and third-harmonic generation, the parallel components are dominant when the emitted harmonic intensity is maximized.
There are qualitative differences between the lowerorder harmonics (the second-and third-order harmonics) and the higher-order harmonics (the fourth-and fifthorder harmonics).In Fig. 2 (c), the fourth harmonic yield reaches a maximum at an angle θ of 90 • under the strongest applied static field, E dc = 1.0 MV/cm.A comparison of Fig. 2 (g) and (k) shows that the perpendicular component dominates the emitted harmonic intensity in this case.In Figs. 2 (d), (h), and (l), the emitted fifth harmonic at the most efficient angle is dominated by the perpendicular component, in spite of the fact that the parallel component is dominant at all angles in the absence of a static field.Hence, the emission paths of the perpendicular components are expected to be important for the enhancement of MIR-induced HHG by a THz field.We observe the same trend for higher-order harmonics (see Appendix B).

C. Comparison of the nonequilibrium steady state and the thermodynamic model
Here, we investigate the role of a nonequilibrium steady state in HHG by comparing the results of the quasistatic approximation and the thermodynamic model [48].As introduced in Sec.III A, the quasistatic approximation consists of replacing the THz pulse with the corresponding static field to describe the electronic system under a THz field.By contrast, the thermodynamic model consists of approximating the electronic system under a THz pulse by a state with a high electron temperature [48].The difference between the quasistatic approximation and the thermodynamic model reflects the difference between the nonequilibrium and equilibrium distributions, clarifying the role of the nonequilibrium steady-state in HHG.
The quasistatic approximation is characterized by the strength of the static field, E dc , whereas the thermodynamic model is characterized by the electron temperature T e .To compare these two models that are formulated using different parameters, we introduce the excess energy [33] as a common measure of the excitation intensity.The excess energy corresponding to the quasistatic approximation is computed as the change in the total energy given in Eq. ( 6) of the nonequilibrium steady state caused by the application of the static field E dc (t).Hence, the excess energy of the nonequilibrium steady-state depends on the static field strength as ∆E non−eq excess (E dc ).The excess energy of the thermodynamic model is computed as the change in the total energy caused by increasing the temperature from room temperature (T e = 300 K) and hence, depends on the electron temperature as ∆E thermo excess (T e ).Thus, converting both the static field strength E dc in the quasistatic approximation and the electron temperature T e in the thermodynamic model to common excess energy enables the two models to be objectively and quantitatively compared [33].
Figure 3 shows the comparison of the results obtained using quasistatic approximation and the thermodynamic model.The MIR field strength is set to 6.5 MV/cm, and the MIR field polarization direction is set to the Γ-M direction (the x-axis of the present setup).The thermodynamic model preserves the intrinsic inversion symmetry of graphene and therefore does not produce even-order harmonics.Hence, we only analyze the odd-order harmonics generated by this model.Figure 3 (a) shows that under a static field, perpendicular and parallel to the MIR polarization, the MIRinduced third harmonic is considerably enhanced and suppressed within the quasistatic approximation but remains almost constant within the thermodynamic model.Figures 3 (b) and (c), show that under a static field, the fifth-and seventh-harmonic yields are significantly enhanced within the quasistatic approximation but decrease slightly as the electron temperature increases within the thermodynamic model.Hence, the enhancement of HHG cannot be described by simple heating of electronic systems within the thermodynamic model and originates from the non-equilibrium nature of fieldinduced electronic dynamics.The small change in the harmonic yields within the thermodynamic model relative to that predicted by the nonequilibrium steady state picture indicates that modification of the population distribution around the Fermi level has little effect on the spectra of HHG.

D. Contribution of the nonequilibrium population
Having demonstrated the importance of the nonequilibrium steady-state under a THz field, we elucidate the role of a coherent coupling between the MIR and THz fields beyond the simple population contribution induced by the THz field.To highlight the coherent coupling contribution, we evaluate the contribution from incoherent coupling by introducing a nonequilibrium population distribution model as an extension of the thermodynamic model.Within the thermodynamic model, the contribution from the THz field is described by modifying the population distribution by increasing the electronic temperature of the reference Fermi-Dirac distribution.Hence, the thermodynamic population only captures the population contribution (the diagonal element of the density matrix) of the THz-induced effect based on the thermal distribution.Here, we extend the thermodynamic model by replacing the reference Fermi-Dirac distribution in the relaxation operator in Eq. ( 4) with the population distribution of the nonequilibrium steady state under a static field.The extended model includes the population contribution (given by diagonal elements of the density matrix) but not THz-induced coherence (given by the off-diagonal elements of the density matrix).Hence, a comparison of the nonequilibrium population model and the fully dynamical model can reveal the contribution from the coherent coupling between the THz and MIR fields.
To formulate the nonequilibrium population model, we first analyze the population distribution in the nonequilibrium steady state under a static field.The population distribution in the Brillouin zone can be expressed as where K (t) is the accelerated wavevector in accordance with the acceleration theorem, K (t) = k + eA(t).
The population distribution in the nonequilibrium steady state can be evaluated in the long-time propagation limit FIG. 3. The emitted light intensity, I nth , is shown as a function of the excess energy for (a) third (b) fifth, and (c) seventh harmonics.The results for the nonequilibrium steady-states induced by a static field parallel (red solid line) and perpendicular (blue dashed line) to the MIR field are compared with the thermodynamic model (green dotted line).In each panel, the field strength of the static field parallel to the MIR field is shown as the secondary axis.
under a static field A(t) = E dc Θ(t), Figure 4 (a) shows the population distribution in the conduction band for the nonequilibrium steady-state under a static field with a strength of E dc = 0.5 MV/cm.The static field is polarized along the Γ-M direction (x-axis).The Dirac point (K point) is depicted by the blue circle.In Fig. 4 (a), the region to the left of the Dirac point is largely occupied by the field-induced population in the nonequilibrium steady-state, whereas the region to the right of the Dirac point is almost empty, breaking the inversion symmetry of the system.We employ this population distribution as the reference distribution of the relaxation operator in Eq. ( 4) instead of the Fermi-Dirac distribution used to construct the nonequilibrium population model in the aforementioned discussion.shows the angular dependence of the second-harmonic yield in the presence of a static field with a strength of E dc = 0.5 MV/cm.The corresponding angular dependences of the third, fourth, and fifth harmonics are shown in Figs. 4 (c-e), respectively.In each panel, the result obtained using the nonequilibrium population model is shown by the blue solid line, and the green solid line corresponds to the result obtained using the fully dynamical model, which is identical to the result shown in Fig. 2. In Figs 4 (b) and (d), the evenorder harmonics computed with the nonequilibrium pop-ulation model are negligibly weak compared with those computed using the fully dynamical calculation.This result indicates that under the charge-neutral condition (µ = 0) investigated here, the two-and four-phonon resonances of the MIR field are far from the Fermi level and cannot be modified by the population changes around the Fermi surface, resulting in small contributions to evenorder harmonic generations.By contrast, within the fully dynamic calculation, the THz field can coherently couple with the MIR field via the off-diagonal elements of the density matrix.Thus, the coherent coupling can be realized both around the Fermi level and anywhere in the Brillouin zone, as long as the dipole transition is allowed.Hence, the coherent coupling contribution may enhance the contribution from the resonant quantum pass, inducing stronger even-order harmonic generation.
Figure 4 (c) shows that for perpendicular THz and MIR fields, the third-harmonic yield calculated using the fully dynamical model is 1.57 times stronger than that computed using the nonequilibrium population model when the two fields are perpendicular.This result indicates that the THz field enhances third-harmonic generation for the perpendicular configuration and that both coherent coupling and the incoherent population play important roles in the THz-assisted enhancement mechanism.By contrast, the third-harmonic yield computed using the fully dynamical calculation is 0.57 times smaller than that computed using the nonequilibrium population model when the two fields are parallel.This result indicates that the contributions from the coherent coupling and the incoherent population cancel each other, weakening the total signal.Therefore, both coherent coupling and incoherent population affect third-harmonic generation under the investigated condition but play different roles depending on the relative angle θ between the THz and MIR fields.
Figure 4 (e) shows that the fifth-order harmonic yield computed using the fully dynamical calculation is considerably larger than that computed using the nonequilibrium population model, except in the range where the THz and MIR fields are parallel.Hence, the coherent coupling is the dominant contribution to the enhancement of fifth-harmonic generation for most angles but the effects of coherent coupling and the incoherent population are both important when the MIR and THz fields are parallel.The consistent results are obtained for higherorder harmonics (see Appendix B).

IV. SUMMARY
We used a quantum master equation to model MIRinduced HHG in graphene in the presence of a strong THz field.We first computed the electron dynamics in graphene by explicitly employing MIR and THz pulses and evaluated the emitted high-harmonic spectra.Next, we developed a quasistatic approximation by analyzing MIR-induced HHG under a static field to replace the THz pulses.The THz-assisted MIR-induced HHG spectra were accurately reproduced by a static field within the quasistatic approximation, thus validating the application of this approximation for describing the induced dynamics generated by an applied strong THz field.
We then investigated the intensity of the emitted harmonics for different relative angles between the MIR and THz fields.In the absence of a THz field, the emitted odd-order harmonics exhibit an almost circular angular dependence, reflecting the circular symmetry of Dirac cones, whereas no even-order harmonics are emitted because of the intrinsic inversion symmetry of graphene.Under an intense THz field, the emitted harmonics exhibit a strong angular dependence along with enhancement and suppression of the harmonic yield.For example, the emitted fifth harmonic can be enhanced 10 times under a THz field with a strength of 0.5 MV/cm with respect to the result without the THz field, as shown in Fig. 2 (d).
To elucidate the mechanism by which a THz field enhances MIR-induced HHG, we compared the results obtained using the quasistatic approximation and the thermodynamic model, which treats the effect of the THz field as a simple increase in the electron temperature of the Fermi-Dirac model [48].The thermodynamic model does not reproduce the enhancement of MIR-induced HHG, indicating that nonequilibrium THz-induced dynamics play an essential role in the enhancement.
To gain further insight into the enhancement of MIRinduced HHG by a THz field, we developed a nonequilibrium population distribution model.Within this model, THz-induced effects are treated as a change in the population distribution in the nonequilibrium steady state.The results obtained using the fully dynamical calculation and the nonequilibrium population distribution model were compared to elucidate the roles of coherent coupling between the MIR and THz fields.The THzinduced even-order harmonics and the THz-enhanced high-order harmonics are dominated by the coherent coupling contribution, whereas the enhancement of the third harmonics under a THz field is affected by both the coherent coupling and the nonequilibrium population.Furthermore, the enhancement of the higher-order harmonics is dominated by the coherent coupling contribution.These enhancement mechanisms are not rigidly limited by the conditions of the laser parameters investigated in this work but can be induced in rather general conditions.Therefore, it is key to control both the coherent coupling and the population for enhancing HHG from solids by employing multicolor laser fields.cellence 'CUI: Advanced Imaging of Matter'-EXC 2056 -project ID 390715994, SFB-925 "Light induced dynamics and control of correlated quantum systems"project 170620586 of the Deutsche Forschungsgemeinschaft (DFG), and the Max Planck-New York City Center for Non-Equilibrium Quantum Phenomena.The Flatiron Institute is a division of the Simons Foundation.This work used computational resources of the HPC systems at the Max Planck Computing and Data Facility (MPCDF) and the Fujitsu PRIMERGY CX400M1/CX2550M5 (Oakbridge-CX) at the Information Technology Center, The University of Tokyo through the HPCI System Research Project (Project ID:hp220112).
Appendix A: Angular dependence of high-order harmonics without THz fields Here, we evaluate the angular dependence of the emitted harmonic yield I nth without THz fields, analyzing the intrinsic symmetry of graphene.Figures 5 show the computed angular dependence of the emitted harmonic yield I nth obtained using only the MIR field in the same conditions as Fig. 2. Reflecting the six-fold symmetry of the hexagonal lattice of graphene, the emitted harmonics also show the six-fold symmetry in the angular dependence.As seen from Fig. 5, the lower-order harmonics exhibit an almost circular angular dependence, reflecting the circular symmetry of Dirac cones.By contrast, the higher-order harmonics exhibit more complex six-fold symmetry in the angular dependence since the electronic structure of graphene deviates from a simple Dirac cone when a single-particle energy is far from the Dirac point.Appendix B: Angular dependence of high-order harmonics Here, we analyze the angular dependence of the highorder harmonics in the same way as that used to analyze Fig. 2. Figures 6 (a) and (b) show the angular dependence of the sixth and seventh harmonics, respectively.Figures 6 (c) and (e) show the sixth-harmonic signal decomposed into parallel and perpendicular components, respectively.The same decomposition is shown for the seventh-order harmonic in Figs. 6 (d) and (f).Consistent with the results for the fourth and fifth harmonics shown in Fig. 2, the perpendicular components make a large contribution to the enhancement of MIRinduced HHG by a THz field, as shown in Fig. 6.Furthermore, we compare the results for the sixth and seventh harmonics obtained using the nonequilibrium population model and the nonequilibrium steady state.Figures 7 (a Consistent with the analysis results shown in Fig. 4, the coherent coupling between the MIR and THz fields plays an essential role in the enhancement of the HHG and goes beyond the simple field-induced population contribution. ) and (b) show the computed current J (t) induced by E THz (t) + E MIR (t) as a function of time.The result with the parallel configuration (e MIR = e x = e THz ) is shown in Fig. 1 (a), while the result with the perpendicular configuration (e MIR = e y ⊥ e THz ) is shown in Fig. 1 (b).As seen from Figs. 1 (a) and (b), the THz field induces a current on the picosecond time scale, whereas the MIR field induces a current on a much shorter time scale.

FIG. 1 .
FIG. 1. (a, b) The current J (t) induced by THz and MIR fields, ETHz(t) + EMIR(t).The inset is the panel (a) shows the time profile of the applied THz field.(c, d) The current J (t) induced by the static and MIR fields, E dc (t) + EMIR(t − τMIR).In the panels (a) and (c), the polarization of all the fields is parallel to the Γ-M direction (the x-direction in the present setup) as eTHz = e dc = eMIR = ex.In the panels (b) and (d), the polarization of THz and static fields is parallel to the x-direction as eTHz = e dc = eMIR = ex, while that of the MIR field is perpendicular as eMIR = ey.(e) The power spectra IHHG(ω) computed using the current in (a) and (c).(f) The power spectra IHHG(ω) computed using the current in (b) and (d).

FIG. 2 .
FIG.2.The angular dependence of the harmonic yield in the nonequilibrium steady-states under a static field along the Γ-M direction is shown.The angle θ denotes the relative angle between the static field and the MIR field.(a-d) The total intensity I nth total is shown for the second, third, fourth, and fifth harmonics.(e-h) The component of the intensity parallel to eMIR is shown for each harmonic.(i-l) The component of the intensity perpendicular to eMIRis shown for each harmonic.The results are normalized by the maximum total intensity I nth total for each harmonic.

FIG. 4 .
FIG. 4. (a) The calculated conduction population distribution, n neq−steady ck for the nonequilibrium steady-state is shown.Here, the Dirac point is indicated by the blue circle.(b-e) The angular dependence of the emitted harmonic intensity is shown for the (b) second, (c) third, (d) fourth, and (e) fifth harmonics.The results obtained using the nonequilibrium population model and the nonequilibrium steady-state are shown by the blue and green solid lines, respectively.

Figure 4 (
Figure 4 (b)  shows the angular dependence of the second-harmonic yield in the presence of a static field with a strength of E dc = 0.5 MV/cm.The corresponding angular dependences of the third, fourth, and fifth harmonics are shown in Figs.4 (c-e), respectively.In each panel, the result obtained using the nonequilibrium population model is shown by the blue solid line, and the green solid line corresponds to the result obtained using the fully dynamical model, which is identical to the result shown in Fig.2.In Figs4 (b) and (d), the evenorder harmonics computed with the nonequilibrium pop-

FIG. 5 .
FIG.5.The angular dependence of the harmonic yield obtained from the electron dynamics calculations in the presence of the MIR field.The numerical conditions are the same as those of the calculations in Fig.2.The third, fifth, and seventh harmonic yields are scaled by factors of 60, 800, and 1000, respectively.

FIG. 6 .
FIG. 6.The angular dependence of the harmonic yields in the nonequilibrium steady state under a static field along the Γ-M direction is shown.The angle θ denotes the relative angle between the static field and the MIR field.(a and b) The total intensity I nth total for the sixth and seventh harmonics is shown, respectively.(c and d) The component of the intensity parallel to eMIR is shown for each harmonic.(e and h) Th component of the intensity perpendicular to eMIR is shown for each harmonic.The results are normalized by the maximum total intensity I nth total of each harmonic.

FIG. 7 .
FIG. 7. The angular dependence of the emitted harmonic intensity for the (a) sixth and (b) seventh harmonics are shown.The results obtained using the nonequilibrium population model and the nonequilibrium steady-state are shown by the blue and green solid lines, respectively.