High-order harmonic generation (HHG) from atomic and molecular gases^[1–4] has been a significant research topic in recent decades due to broad applications, such as diffractive imaging,^[5] extracting of structural information of atoms or molecules,^[6–12] probing magnetic dynamics.^[13] More importantly, we can construct untrashort attosecond (as) pulses from harmonics,^[14–19] which is an important tool for observing and controlling the ultrafast electronic dynamics in atoms and molecules. The physical picture of the atomic HHG process can be easily understood by the semiclassical three-step model.^[20] According to this model, the cutoff energy of harmonic spectrum is I_p + 3.17U_p, where I_p is the atomic ionization potential, U_p is ponderomotive energy which is proportional to the intensity of laser pulse. Very recently, HHG from the crystalline solid^[21–28] has also aroused great interests of researchers.^[29–57] In this case, both experimental and theoretical investigations show that the harmonic cutoff energy of the first plateau increases linearly with the amplitude of field field. To interpret the mechanism of harmonic emission from crystals, a similar three-step model in momentum space is proposed on a basis of the energy band structure,^[30] which can describe well the motion of electron-hole pairs in crystals: firstly, the electrons in the valence band driven by laser field tunnel to the conduction band, producing the electron-hole pairs; then, the electron-hole pairs oscillate in crystal lattice under laser filed and gain momentum (energy); finally, the electrons recombine with holes and radiate photons (HHG) later. Obviously, for crystal matters, the strength of HHG depends on the tunneling probability between valence band and conduction band, which is related to the energy gap and strength of laser, according to the Landau–Zaner law.

Theoretically, the semiconductor Bloch equations (SBEs) are often adopted to investigate the HHG emission from crystal, for it can easily introduce dephase time phenomenologically to take the relaxation effect into account.^[35] However, for crystals with the zero band gap, such as graphene, when the initial electron momentum equals 0, the dipole moment will be divergent near the Dirac points in the length gauge. Thus, one must dig out the Dirac point or take a constant in the calculations,^[58,59] which is inconvenient and might influence the accuracy of calculations. Alternatively, a simple transformation^[60,61] is used to deal with the divergence due to the singular values of dipole moment near Dirac points. Furthermore, we could also solve directly time-dependent Dirac equation (TDDE).^[62] However, in the last two cases it is difficult to consider the dephase time.

In this work, we present a method of solving SBEs in velocity gauge, which can avoid the divergence of dipole transition moment near Dirac points in a single-layer graphene. To confirm our method, we first compare our predictions for the population and electron current with other methods. Using our method, next we analyze the contributions of quantum trajectories around Dirac cone, and investigate the behavior of harmonic from graphene generated by elliptically polarized laser pulses. Finally we study HHG dependence on crystallographic orientation of graphene. The article is organized as follows. In Section 2, we give our theoretical framework. In Section 3, we show numerical simulations on HHG. Then we conclude our work in Section 4.

Our theoretical approach is based on the time-dependent tight-binding (TB) approximation,^[63–65] which has been extensively applied in the investigation of laser–graphene interaction.^[58,59,66] For the graphene, the most important electronic states responsible for the harmonic emission are those arising from the carbon 2p_z orbitals, which form the extended π-bands. With the consideration of nearest-neighbor tight-binding of the π-electrons, we can easily obtain the energy bands. Near the Dirac points, energy bands calculated with T–B approximation are in accordance well with those obtained from first-principle calculations, i.e., density functional theory, as long as the laser intensity is not too higher. The time-dependent Hamiltonian of graphene in the length gauge can be written as^[67]

$$ \begin{eqnarray}\hat{H}(t)=\left[\begin{array}{cc}0 & \gamma f({\boldsymbol{k}}(t))\\ \gamma {f}^{* }({\boldsymbol{k}}(t)) & 0\end{array}\right],\end{eqnarray}$$ (1)

In order to avoid this divergence, we can calculate current in velocity gauge with Bloch basis. We can write H^(t) as

$$ \begin{eqnarray}\hat{H}(t)={\hat{H}}_{0}+{\boldsymbol{A}}(t)\cdot {\boldsymbol{P}}+{\boldsymbol{A}}{(t)}^{2}=\displaystyle \sum _{\alpha =1}^{3}{\hat{H}}_{0}^{\alpha }+\displaystyle \sum _{\alpha =1}^{3}{\hat{H}}_{{\rm{I}}}^{\alpha },\end{eqnarray}$$ (6)

$$ \begin{eqnarray}\begin{array}{lll}{\hat{H}}_{0} & = & \gamma \left[\begin{array}{cc}0 & {{\rm{e}}}^{{\rm{i}}{\boldsymbol{k}}{\delta }_{1}}\\ {{\rm{e}}}^{-{\rm{i}}{\boldsymbol{k}}{\delta }_{1}} & 0\end{array}\right]+\gamma \left[\begin{array}{cc}0 & {{\rm{e}}}^{{\rm{i}}{\boldsymbol{k}}{\delta }_{2}}\\ {{\rm{e}}}^{-{\rm{i}}{\boldsymbol{k}}{\delta }_{2}} & 0\end{array}\right]\\ & & +\,\gamma \left[\begin{array}{cc}0 & {{\rm{e}}}^{{\rm{i}}{\boldsymbol{k}}{\delta }_{3}}\\ {{\rm{e}}}^{-{\rm{i}}{\boldsymbol{k}}{\delta }_{3}} & 0\end{array}\right],\end{array}\end{eqnarray}$$ (7)

$$ \begin{eqnarray}\begin{array}{lll}{\hat{H}}_{I}^{\alpha }(t)& = & \gamma \left[\begin{array}{cc}0 & {{\rm{e}}}^{{\rm{i}}{\boldsymbol{k}}{\delta }_{\alpha }}\\ {{\rm{e}}}^{-{\rm{i}}{\boldsymbol{k}}{\delta }_{\alpha }} & 0\end{array}\right]\\ & & \times \left[\begin{array}{cc}{{\rm{e}}}^{-{\rm{i}}{\boldsymbol{A}}(t){\delta }_{\alpha }}-1 & 0\\ 0 & {{\rm{e}}}^{{\rm{i}}{\boldsymbol{A}}(t){\delta }_{\alpha }}-1\end{array}\right].\end{array}\end{eqnarray}$$ (8)

Utilizing Bloch basis as unitary transformation matrix

$$ \begin{eqnarray}U=\displaystyle \frac{1}{\sqrt{2}}\left[\begin{array}{cc}{{\rm{e}}}^{{\rm{i}}{\theta }_{f({\boldsymbol{k}})}} & -{{\rm{e}}}^{{\rm{i}}{\theta }_{f({\boldsymbol{k}})}}\\ 1 & 1\end{array}\right],\end{eqnarray}$$ (9)

 1

$$ \begin{eqnarray}\begin{array}{lll} & & U{\hat{H}}_{0}^{\alpha }{U}^{\dagger }=U \left[\begin{array}{cc}0 & \gamma {{\rm{e}}}^{{\rm{i}}{\boldsymbol{k}}{\delta }_{\alpha }}\\ \gamma {{\rm{e}}}^{-{\rm{i}}{\boldsymbol{k}}{\delta }_{\alpha }} & 0\end{array}\right]{U}^{\dagger },\end{array}\end{eqnarray}$$ (10)

$$ \begin{eqnarray}\begin{array}{lll} & & U{\hat{H}}_{I}^{\alpha }{U}^{\dagger }=U \left[\begin{array}{cc}0 & \gamma {{\rm{e}}}^{{\rm{i}}{\boldsymbol{k}}{\delta }_{\alpha }}\\ \gamma {{\rm{e}}}^{-{\rm{i}}{\boldsymbol{k}}{\delta }_{\alpha }} & 0\end{array}\right]{U}^{\dagger }\\ & & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \times U \left[\begin{array}{cc}{{\rm{e}}}^{-{\rm{i}}{\boldsymbol{A}}(t){\delta }_{\alpha }}-1 & 0\\ 0 & {{\rm{e}}}^{{\rm{i}}{\boldsymbol{A}}(t){\delta }_{\alpha }}-1\end{array}\right]{U}^{\dagger }.\ \ \ \ \ \ \ \ \end{array}\end{eqnarray}$$ (11)

$$ \begin{eqnarray}{\hat{H}}_{B}(t)=U\hat{H}(t){U}^{\dagger }=\gamma \left[\begin{array}{cc}{W}_{1}(t) & {W}_{2}(t)\\ {W}_{3}(t) & {W}_{4}(t)\end{array}\right]\end{eqnarray}$$ (12)

$$ \begin{eqnarray}\begin{array}{lll} & & {W}_{1}(t)=-f({\boldsymbol{k}}(t)){{\rm{e}}}^{-{\rm{i}}{\theta }_{f({\boldsymbol{k}})}}-{f}^{* }({\boldsymbol{k}}(t)){{\rm{e}}}^{{\rm{i}}{\theta }_{f({\boldsymbol{k}})}},\end{array}\end{eqnarray}$$ (13)

$$ \begin{eqnarray}\begin{array}{lll} & & {W}_{2}(t)=-f({\boldsymbol{k}}(t)){{\rm{e}}}^{-{\rm{i}}{\theta }_{f({\boldsymbol{k}})}}+{f}^{* }({\boldsymbol{k}}(t)){{\rm{e}}}^{{\rm{i}}{\theta }_{f({\boldsymbol{k}})}},\end{array}\end{eqnarray}$$ (14)

$$ \begin{eqnarray}\begin{array}{lll} & & {W}_{3}(t)=f({\boldsymbol{k}}(t)){{\rm{e}}}^{-{\rm{i}}{\theta }_{f({\boldsymbol{k}})}}-{f}^{* }({\boldsymbol{k}}(t)){{\rm{e}}}^{{\rm{i}}{\theta }_{f({\boldsymbol{k}})}},\end{array}\end{eqnarray}$$ (15)

$$ \begin{eqnarray}\begin{array}{lll} & & {W}_{4}(t)=f({\boldsymbol{k}}(t)){{\rm{e}}}^{-{\rm{i}}{\theta }_{f({\boldsymbol{k}})}}+{f}^{* }({\boldsymbol{k}}(t)){{\rm{e}}}^{{\rm{i}}{\theta }_{f({\boldsymbol{k}})}},\end{array}\end{eqnarray}$$ (16)

$$ \begin{eqnarray}\begin{array}{lll} & & \displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{C}_{{\rm{v}}}=\displaystyle \frac{\gamma }{2}{\rm{i}}[-{W}_{1}(t){C}_{{\rm{v}}}-{W}_{2}(t){C}_{{\rm{c}}}],\end{array}\end{eqnarray}$$ (17)

$$ \begin{eqnarray}\begin{array}{lll} & & \displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{C}_{{\rm{c}}}=\displaystyle \frac{\gamma }{2}{\rm{i}}[-{W}_{3}(t){C}_{{\rm{v}}}-{W}_{4}(t){C}_{{\rm{c}}}],\end{array}\end{eqnarray}$$ (18)

Next following the Ref. [55], we obtain following dynamic equations

$$ \begin{eqnarray}\begin{array}{lll} & & \displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\rho ({\boldsymbol{k}},t)=-{\rm{i}}{A}_{4}^{* }\rho ({\boldsymbol{k}},t)+{\rm{i}}{A}_{1}^{* }\rho ({\boldsymbol{k}},t)\\ & & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-{\rm{i}}{A}_{3}(1-{f}_{{\rm{e}}}-{f}_{{\rm{h}}})-{\gamma }_{{\rm{r}}}\rho ({\boldsymbol{k}},t),\end{array}\end{eqnarray}$$ (19)

$$ \begin{eqnarray}\begin{array}{lll} & & \displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{f}_{{\rm{e}}}=2\,{\rm{Im}}\,[{A}_{3}^{* }\rho ({\boldsymbol{k}},t)]-{\gamma }_{{\rm{l}}}{f}_{{\rm{e}}},\end{array}\end{eqnarray}$$ (20)

$$ \begin{eqnarray}\begin{array}{lll} & & \displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{f}_{{\rm{h}}}=2\,{\rm{Im}}\,[-{A}_{2}^{* }\rho ({\boldsymbol{k}},t)]-{\gamma }_{{\rm{l}}}{f}_{{\rm{h}}},\end{array}\end{eqnarray}$$ (21)

$$ \begin{eqnarray}\begin{array}{lll} & & {A}_{1}(t)=-\displaystyle \frac{\gamma }{2}[f({\boldsymbol{k}}(t)){{\rm{e}}}^{-{\rm{i}}{\theta }_{f({\boldsymbol{k}})}}+{f}^{* }({\boldsymbol{k}}(t)){{\rm{e}}}^{{\rm{i}}{\theta }_{f({\boldsymbol{k}})}}],\end{array}\end{eqnarray}$$ (22)

$$ \begin{eqnarray}\begin{array}{lll} & & {A}_{2}(t)=-\displaystyle \frac{\gamma }{2}[f({\boldsymbol{k}}(t)){{\rm{e}}}^{-{\rm{i}}{\theta }_{f({\boldsymbol{k}})}}-{f}^{* }({\boldsymbol{k}}(t)){{\rm{e}}}^{{\rm{i}}{\theta }_{f({\boldsymbol{k}})}}],\end{array}\end{eqnarray}$$ (23)

$$ \begin{eqnarray}\begin{array}{lll} & & {A}_{3}(t)=\displaystyle \frac{\gamma }{2}[f({\boldsymbol{k}}(t)){{\rm{e}}}^{-{\rm{i}}{\theta }_{f({\boldsymbol{k}})}}-{f}^{* }({\boldsymbol{k}}(t)){{\rm{e}}}^{{\rm{i}}{\theta }_{f({\boldsymbol{k}})}}],\end{array}\end{eqnarray}$$ (24)

$$ \begin{eqnarray}\begin{array}{lll} & & {A}_{4}(t)=\displaystyle \frac{\gamma }{2}[f({\boldsymbol{k}}(t)){{\rm{e}}}^{-{\rm{i}}{\theta }_{f({\boldsymbol{k}})}}+{f}^{* }({\boldsymbol{k}}(t)){{\rm{e}}}^{{\rm{i}}{\theta }_{f({\boldsymbol{k}})}}].\end{array}\end{eqnarray}$$ (25)

$$ \begin{eqnarray}\begin{array}{lll} & & \displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{\rho }_{{\boldsymbol{k}}}=-{\rm{i}}{B}_{4}{\rho }_{{\boldsymbol{k}}}+{\rm{i}}{B}_{1}{\rho }_{{\boldsymbol{k}}}-{\rm{i}}{B}_{3}(1-{f}_{{\rm{e}}}-{f}_{{\rm{h}}})-{\gamma }_{{\rm{r}}}{\rho }_{{\boldsymbol{k}}},\ \ \ \ \ \ \ \ \end{array}\end{eqnarray}$$ (26)

$$ \begin{eqnarray}\begin{array}{lll} & & \displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{f}_{{\rm{e}}}=2\,{\rm{iIm}}\,({B}_{3}^{* }{\rho }_{{\boldsymbol{k}}})-{\gamma }_{{\rm{l}}}{f}_{{\rm{e}}},\end{array}\end{eqnarray}$$ (27)

$$ \begin{eqnarray}\begin{array}{lll} & & \displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{f}_{{\rm{h}}}=-2\,{\rm{iIm}}\,({B}_{2}^{* }{\rho }_{{\boldsymbol{k}}})-{\gamma }_{{\rm{l}}}{f}_{{\rm{h}}},\end{array}\end{eqnarray}$$ (28)

$$ \begin{eqnarray}\begin{array}{lll} & & {B}_{1}(t)={E}_{v}({\boldsymbol{k}})-{v}_{{\rm{f}}}\cos ({\theta }_{{\boldsymbol{k}}}){A}_{x}-{v}_{{\rm{f}}}\sin ({\theta }_{{\boldsymbol{k}}}){A}_{y},\end{array}\end{eqnarray}$$ (29)

$$ \begin{eqnarray}\begin{array}{lll} & & {B}_{2}(t)=-{\rm{i}}{v}_{{\rm{f}}}\sin ({\theta }_{{\boldsymbol{k}}}){A}_{x}+{\rm{i}}{v}_{{\rm{f}}}{A}_{y}\cos ({\theta }_{{\boldsymbol{k}}}),\end{array}\end{eqnarray}$$ (30)

$$ \begin{eqnarray}\begin{array}{lll} & & {B}_{3}(t)={\rm{i}}{v}_{{\rm{f}}}\sin ({\theta }_{{\boldsymbol{k}}}){A}_{x}-{\rm{i}}{v}_{{\rm{f}}}{A}_{y}\cos ({\theta }_{{\boldsymbol{k}}}),\end{array}\end{eqnarray}$$ (31)

$$ \begin{eqnarray}\begin{array}{lll} & & {B}_{4}(t)={E}_{{\rm{c}}}({\boldsymbol{k}})+{v}_{{\rm{f}}}\cos ({\theta }_{{\boldsymbol{k}}}){A}_{x}+{v}_{{\rm{f}}}\sin ({\theta }_{{\boldsymbol{k}}}){A}_{y},\end{array}\end{eqnarray}$$ (32)

The above equation sets can be numerically solved for each independent k by the classical fourth-order Runge–Kutta method combining with an adaptive step-size routine. We then obtain the single-electron current. Around K, we have

$$ \begin{eqnarray}\begin{array}{lll} & & {j}_{x}^{K}(t)=(1-{f}_{{\rm{e}}}-{f}_{{\rm{h}}})\cos ({\theta }_{{\boldsymbol{k}}})-2{\rm{Im}}({\rho }_{{\boldsymbol{k}}})\sin ({\theta }_{{\boldsymbol{k}}}),\end{array}\end{eqnarray}$$ (33)

$$ \begin{eqnarray}\begin{array}{lll} & & {j}_{y}^{K}(t)=(1-{f}_{{\rm{e}}}-{f}_{{\rm{h}}})\sin ({\theta }_{{\boldsymbol{k}}})+2{\rm{Im}}({\rho }_{{\boldsymbol{k}}})\cos ({\theta }_{{\boldsymbol{k}}}).\ \ \ \ \ \ \end{array}\end{eqnarray}$$ (34)

$$ \begin{eqnarray}\begin{array}{lll} & & {j}_{x}^{{K}^{\prime}}(t)=(1-{f}_{{\rm{e}}}-{f}_{{\rm{h}}})\cos ({\theta }_{{\boldsymbol{k}}})+2{\rm{Im}}({\rho }_{{\boldsymbol{k}}})\sin ({\theta }_{{\boldsymbol{k}}}),\end{array}\end{eqnarray}$$ (35)

$$ \begin{eqnarray}\begin{array}{lll} & & {j}_{y}^{{K}^{\prime}}(t)=(1-{f}_{{\rm{e}}}-{f}_{{\rm{h}}})\sin ({\theta }_{{\boldsymbol{k}}})-2{\rm{Im}}({\rho }_{{\boldsymbol{k}}})\cos ({\theta }_{{\boldsymbol{k}}}).\ \ \ \ \ \ \end{array}\end{eqnarray}$$ (36)

After we taking into account the Fermi distribution and initial condition 1 − f_e − f_h = −F(p) + F(−p), where F(p) = 1 + exp[(E_n(p) − u)/k_BT]⁻¹, with u, k_B, and T being the chemical potential, Boltzman constant, and temperature, respectively, we obtain the total integrated electric current J(t). The Harmonic spectrum is given by I(ω) ∝ |ω J(ω)|², J(ω) is Fourier transform of the J(t).

If we want to further gain the information from valence and conduction bands, at each time, we can easily get density-matrix elements ρ_mn by simple unitary transformation

$$ \begin{eqnarray}\begin{array}{lll}{\rho }_{cc} & = & 0.25({Z}_{3}^{* }{Z}_{3})(1-{f}_{{\rm{e}}})+0.25({Z}_{3}^{* }{Z}_{4}){\rho }_{{\boldsymbol{k}}}\\ & & +\,0.25({Z}_{4}^{* }{Z}_{3}){\rho }_{{\boldsymbol{k}}}^{* }+0.25({Z}_{4}^{* }{Z}_{4}){f}_{{\rm{h}}},\end{array}\end{eqnarray}$$ (37)

$$ \begin{eqnarray}\begin{array}{lll}{\rho }_{cv} & = & 0.25({Z}_{3}^{* }{Z}_{1})(1-{f}_{{\rm{e}}})+0.25({Z}_{3}^{* }{Z}_{2}){\rho }_{{\boldsymbol{k}}}\\ & & +\,0.25({Z}_{4}^{* }{Z}_{1}){\rho }_{{\boldsymbol{k}}}^{* }+0.25({Z}_{4}^{* }{Z}_{2}){f}_{{\rm{h}}},\end{array}\end{eqnarray}$$ (38)

$$ \begin{eqnarray}\begin{array}{lll}{\rho }_{vv} & = & 0.25({Z}_{1}^{* }{Z}_{1})(1-{f}_{{\rm{e}}})+0.25({Z}_{1}^{* }{Z}_{2}){\rho }_{{\boldsymbol{k}}}\\ & & +\,0.25({Z}_{2}^{* }{Z}_{1}){\rho }_{{\boldsymbol{k}}}^{* }+0.25({Z}_{2}^{* }{Z}_{2}){f}_{{\rm{h}}},\end{array}\end{eqnarray}$$ (39)

$$ \begin{eqnarray}\begin{array}{lll} & & {Z}_{1}(t)=({{\rm{e}}}^{-{\rm{i}}{\theta }_{f({\boldsymbol{k}}(t))}}{{\rm{e}}}^{-{\rm{i}}{\theta }_{f({\boldsymbol{k}})}}+1)/2,\end{array}\end{eqnarray}$$ (40)

$$ \begin{eqnarray}\begin{array}{lll} & & {Z}_{2}(t)=(-{{\rm{e}}}^{-{\rm{i}}{\theta }_{f({\boldsymbol{k}}(t))}}{{\rm{e}}}^{-{\rm{i}}{\theta }_{f({\boldsymbol{k}})}}+1)/2,\end{array}\end{eqnarray}$$ (41)

$$ \begin{eqnarray}\begin{array}{lll} & & {Z}_{3}(t)={Z}_{2}(t),\end{array}\end{eqnarray}$$ (42)

$$ \begin{eqnarray}\begin{array}{lll} & & {Z}_{4}(t)={Z}_{1}(t).\end{array}\end{eqnarray}$$ (43)

From the above equations, we see that there are no divergence around Dirac points in our method.

In the following, we aim to confirm the validity of our method. We first compare our calculated conduction band population with those from Ref. [58], in which Kelardeh et al. investigated the charge transfer process in graphene by a two-band model based on Houston basis. The laser pulse has the form of

$$ \begin{eqnarray}F(t)={F}_{0}{{\rm{e}}}^{-{u}^{2}}(1-2{u}^{2}),\end{eqnarray}$$ (44)

One can see that the electrons movement under laser field lead to interference fringes, and at the end of pulse, the distribution becomes completely symmetric which can be modulated by laser amplitude. Our results with SBEs in the velocity gauge agree very well with those from the two-band model, for both amplitudes, i.e., F₀ = 0.8 V/Å and a higher one of 2.25 V/Å.

In Fig. 3, we further compare ρ_cc(k,t) obtained with two methods for t = 0.75 fs and t = 2.25 fs. It shows again that our results are in good accordance with those from two-band model.^[58]

Next, we compare current calculated by our method (Eqs. (33)–(36)) with those from solving time-dependent Dirac equation (TDDE).^[62] The laser field used in the simulation is a flat-top pulse with a half-cycle ramp on

$$ \begin{eqnarray}{\boldsymbol{A}}(t)=\left\{\begin{array}{ll}0, & {\rm{if}}\,({\rm{t}}\lt 0),\\ {A}_{0}\displaystyle \frac{\omega t}{\pi }\sin (\omega t), & {\rm{if}}\,(0\le {\rm{t}}\lt \pi /\omega ),\\ {A}_{0}\sin (\omega t), & {\rm{if}}\,(\pi /\omega \le {\rm{t}}),\end{array}\right.\end{eqnarray}$$ (45)

We note that another advantage of our method is that in our simulations without solving the gradient of d_cv(k),^[55] which promotes calculation efficiency greatly. In addition, if we use the Houston basis, it is difficult to obtain correct current near Dirac points, although it provides an advantage of decoupling the states at different crystal momentum. From Fig. 4(a) we see that a step transition occurred in the j_x calculated with the Houston basis, which is unphysical.

After proving our method is feasible and effective, we show a few examples of its applications in simulation of HHG from graphene. Firstly, we strive to extract the contributions of individual quantum trajectories around single Dirac cone, which is important in analyzing electron dynamic motion. To achieve this we use the windowed Fourier transform (WFT), form of the wave function is given in Ref. [70]. In that paper, Chizhova et al. showed that the individual quantum trajectories can be extracted by solving TDDE.

In Fig. 5, we show the laser vector potential and generated electron current, wave packets for three initial energies. For initial energy of |E_n| < v_fA₀ and zero p_y, one can see that electron moves in valence band or conduction band driven by laser pulse with a constant current [Figs. 5(a) and 5(b)], the response features intraband dynamics.

Considering an electron with initial momentum p_x = 0.8A₀, p_y = 0.1A₀, it moves initially in the valence band. When laser field reaches its maximum after one optical cycle, electron is driven up to the Dirac point, then the electronic wave packet splits into two parts, one is excited to the conduction band via inter-band Landau-Zener tunneling, while the other remains in the valence band [see Fig. 5(d)]. The superposition of two wave packets results in high-frequency oscillations of the current j_x [Fig. 5(c)], which is determined by coherent superposition of inter-band and intra-band dynamics. The strength of Landau–Zener tunneling rate depends exponentially on the energy gap between the two cones at the p_y = constant conical intersection. The highest oscillation frequency is realized when the two paths reach their maximal separation in energy (near 1.5T).

When initial energies |E_n| ≈ v_fA₀ and p_y ≈ A₀, the total current just shows intra-band dynamics, and electrons move in valence band under laser pulses due to a big energy gap.

Unlike atom and uniaxial hexagonal crystal, HHG from graphene shows abnormal dependence of harmonic strength on the polarization of driving laser field.^[71] Now we move to this issue. We simulate harmonic spectra with the following laser pulse:

$$ \begin{eqnarray}{\boldsymbol{F}}(t)=\displaystyle \frac{{F}_{0}}{\sqrt{1+{\varepsilon }^{2}}}f(t)[\cos (\omega t){\hat{e}}_{x}-\varepsilon \sin (\omega t){\hat{e}}_{y}]\end{eqnarray}$$ (46)

In Fig. 7, we show more detailed dependence of intensity of harmonics 3rd, 5th, 7th, and 9th on the ellipticity of driving laser pulse. It is found that this dependence strongly relies on the harmonic radiated direction. In the y direction, the 3rd harmonic emission has a weak dependence on ε and reaches maximal value at ε ≈ 0.15, while strong radiation in the x direction, and the total harmonic yields decrease monotonically with the increasing of parameter ε. This trend is similar to the case of harmonic radiated from atoms.^[72] For the 5th, 7th, and 9th harmonics, their strength in two directions intersect at one laser ellipticity, and the x-direction harmonic intensity decays quickly with the change of parameter ε, while the y-direction harmonics increases to a maximal value at a certain ε (for example, ε ≈ 0.3 for the 5th, 7th harmonics, ε ≈ 0.2 for the 9th harmonics) and then decreases monotonically with further increasing of ε. It is worth mentioning that the total harmonic yields will display a non-monotonic variation, which can well match the recent experimental observation^[71] and theoretical simulations.^[59]

The polarization direction of laser pulse with respect to the crystal axis is another factor that influence the harmonic emission.^[73] However, the orientation dependence of the 7th harmonic does not reflect the hexagonal structure of graphene in an experimental study.^[71] Here we reinvestigate this issue using our method. We simulate harmonic spectra with a linear laser field as described by Eq. (46) with F₀ = 8.7 × 10⁶ V/cm, ε = 0, and the laser duration is 12 optical cycles. Figure 8 shows our main results.

From Fig. 8(a), we can see that for lower harmonics, same as experimental findings,^[71] their strength almost has no dependence on the laser rotation angle, while for the higher harmonics, their intensities change periodically twice with the increasing the θ from 0° to 120°, which indicates sixfold symmetry structure of graphene.

These results can be understood as follows. Under a weaker laser field, electrons moves mainly around Dirac cones, and these electrons will contribute to lower order harmonics, so band structure of a solid appears unimportant. With the increase of field strength, electrons are driven beyond Dirac cones, inter-band current will play major role and give rise to higher harmonics [Figs. 8(c) and 8(d)], and results reveal a strong link between sixfold symmetry structure of graphene and inter-band polarization [Fig. 8(b)].

In summary, in order to simulate the HHG from graphene we propose a new efficient theoretical approach in this paper. Unlike the previous methods, we solve the semiconductor Bloch equations in velocity gauge, which has the advantages of no defect of divergence of dipole transition moment around Dirac points, avoiding the solve the gradient of dipole transition moment.

Using our method, we have investigated some typical issues associated with HHG from graphene. Our simulations demonstrate that our method can be used successfully to analyze the electron dynamics in valence band and conduction band, especially for the initial condition p_y = 0. We show that harmonic yields in two perpendicular directions depend differently on the ellipticity of driving laser pulse, and the 3rd-harmonic intensity decreases monotonically with the laser ellipticity, while the 5th, 7th, and 9th harmonics can be enhanced at a particular ellipticity, and it reaches to its maximum when ε ≈ 0.3, which are in agreement with the experimental findings. Our simulations for the HHG dependence on the crystal orientation also show that sixfold symmetry of graphene can be retrieved from higher harmonics, which can be further confirmed in the experiments.

In addition to graphene, our method can be also used to investigate the interaction of strong laser fields with other zero-band-gap systems, such as surface states of topological insulators.^[74,75]