Higher Harmonic Generation of Coherent Sub-THz Radiation at Lifshitz Transition in Gapped Bilayer Graphene

Abstract

The microscopic quantum theory of nonlinear interaction of strong coherent electromagnetic radiation with gapped bilayer graphene is used for consideration of the high harmonic generation at low-energy photon excitation Lifshitz transition. The Liouville–von Neumann equation for the density matrix is solved numerically at the nonadiabatic multiphoton excitation regime. By numerical solutions, we examine the rates of the second and third harmonics generation at the particle-hole annihilation at Lifshitz transitions in the linearly polarized coherent electromagnetic wave. The obtained results show that the gapped bilayer graphene can serve as an effective medium for the generation of even and odd high harmonics in the sub-THz domain of frequencies.

APPENDIX

The Liouville–von Neumann equation for a single-particle density matrix can be presented by the form

$${{\rho }_{{\alpha ,\beta }}}({\mathbf{p}},t) = \langle \hat {a}_{{{\mathbf{p}},\beta }}^{\dag }(t){{\hat {a}}_{{{\mathbf{p}},\alpha }}}(t)\rangle ,$$

(25)

where ${{\hat {a}}_{{{\mathbf{p}},\alpha }}}$(t) obeys the Heisenberg equation

$$i\hbar \frac{{\partial {{{\hat {a}}}_{{{\mathbf{p}},\alpha }}}(t)}}{{\partial t}} = [{{\hat {a}}_{{{\mathbf{p}},\alpha }}}(t),\hat {H}].$$

(26)

Due to the homogeneity of the problem we only need the p-diagonal elements of the density matrix. So, taking into account Eqs. (11)–(26), the evolutionary-equation will be

$$\begin{gathered} i\hbar \frac{{\partial {{\rho }_{{\alpha ,\beta }}}({\mathbf{p}},t)}}{{\partial t}} - i\hbar e{\mathbf{E}}(t)\frac{{\partial {{\rho }_{{\alpha ,\beta }}}({\mathbf{p}},t)}}{{\partial {\mathbf{p}}}} \\ = ({{\mathcal{E}}_{\alpha }}({\mathbf{p}}) - {{\mathcal{E}}_{\beta }}({\mathbf{p}}) - i\hbar \Gamma (1 - {{\delta }_{{\alpha \beta }}})){{\rho }_{{\alpha ,\beta }}}({\mathbf{p}},t) \\ + \,{\mathbf{E}}(t)({{D}_{m}}(\alpha ,{\mathbf{p}}) - {{D}_{m}}(\beta ,{\mathbf{p}})){{\rho }_{{\alpha ,\beta }}}({\mathbf{p}},t) \\ + \,{\mathbf{E}}(t)[{{D}_{t}}(\alpha ,{\mathbf{p}}){{\rho }_{{ - \alpha ,\beta }}}({\mathbf{p}},t) - {{D}_{t}}( - \beta ,{\mathbf{p}}){{\rho }_{{\alpha , - \beta }}}({\mathbf{p}},t)], \\ \end{gathered} $$

(27)

nondiagonal Γ—the damping rate. In Eq. (27) the nondiagonal elements are interband polarization ρ_{1, –1}(p, t) = P(p, t) and its complex conjugate ρ_{–1, 1}(p, t) = P*(p, t), and the diagonal elements represent particle distribution functions for conduction N_c(p, t) = ρ_1,1(p, t) and valence ${{N}_{{v}}}$(p, t) = ρ_–1,1(p, t) bands. We will solve the set of differential equations for these functions:

$$\begin{gathered} i\hbar \frac{{\partial {{N}_{c}}({\mathbf{p}},t)}}{{\partial t}} - i\hbar e{\mathbf{E}}(t)\frac{{\partial {{N}_{c}}({\mathbf{p}},t)}}{{\partial {\mathbf{p}}}} \\ = {\mathbf{E}}(t){{{\mathbf{D}}}_{t}}({\mathbf{p}})P{\text{*}}({\mathbf{p}},t) - {\mathbf{E}}(t){\mathbf{D}}_{t}^{*}({\mathbf{p}})P({\mathbf{p}},t), \\ \end{gathered} $$

(28)

$$\begin{gathered} i\hbar \frac{{\partial {{N}_{{v}}}({\mathbf{p}},t)}}{{\partial t}} - i\hbar e{\mathbf{E}}(t)\frac{{\partial {{N}_{{v}}}({\mathbf{p}},t)}}{{\partial {\mathbf{p}}}} \\ = - {\mathbf{E}}(t){{{\mathbf{D}}}_{t}}({\mathbf{p}})P{\text{*}}({\mathbf{p}},t) + {\mathbf{E}}(t){\mathbf{D}}_{t}^{*}({\mathbf{p}})P({\mathbf{p}},t), \\ \end{gathered} $$

(29)

$$\begin{gathered} i\hbar \frac{{\partial P({\mathbf{p}},t)}}{{\partial t}} - i\hbar e{\mathbf{E}}(t)\frac{{\partial P({\mathbf{p}},t)}}{{\partial {\mathbf{p}}}} \\ = [2{{\mathcal{E}}_{1}}({\mathbf{p}}) + {\mathbf{E}}(t){{{\mathbf{D}}}_{m}}({\mathbf{p}}) - i\hbar \Gamma ]P({\mathbf{p}},t) \\ + \,{\mathbf{E}}(t){{{\mathbf{D}}}_{t}}({\mathbf{p}})[{{N}_{{v}}}({\mathbf{p}},t) - {{N}_{c}}({\mathbf{p}},t)]. \\ \end{gathered} $$

(30)

The total mean dipole moments are

$$\begin{gathered} {{D}_{{xm}}}({\mathbf{p}}) = - \frac{{e\hbar U}}{{2{{\mathcal{E}}_{1}}({\mathbf{p}})(\mathcal{E}_{1}^{2}({\mathbf{p}}) - {{U}^{2}}{\text{/}}4)}} \\ \times \left[ {\left( {\frac{{{{p}^{2}}}}{{2m}} - m{v}_{3}^{2}} \right)\frac{{\zeta {{p}_{y}}}}{m} + \frac{{{{{v}}_{3}}}}{m}{{p}_{x}}{{p}_{y}}} \right], \\ \end{gathered} $$

(31)

$$\begin{gathered} {{D}_{{ym}}}({\mathbf{p}}) - \frac{{e\hbar U}}{{2{{\mathcal{E}}_{1}}({\mathbf{p}})(\mathcal{E}_{1}^{2}({\mathbf{p}}) - {{U}^{2}}{\text{/}}4)}} \\ \times \left[ {\left( { - \frac{{{{p}^{2}}}}{{2m}} + m{v}_{3}^{2}} \right)\frac{{\zeta {{p}_{x}}}}{m} + \frac{{{{{v}}_{3}}}}{{2m}}(p_{x}^{2} - p_{y}^{2})} \right]. \\ \end{gathered} $$

(32)

The components of the transition dipole moments are calculated via Eq. (14) by spinor wave functions (6) (see also in [30, 33]):

$$\begin{gathered} {{D}_{{tx}}}({\mathbf{p}}) = \frac{{e\hbar }}{{2{{\mathcal{E}}_{1}}({\mathbf{p}})\sqrt {\mathcal{E}_{1}^{2}({\mathbf{p}}) - {{U}^{2}}{\text{/}}4} }} \\ \times \left( {\left[ {\left( {\frac{{{{p}^{2}}}}{{2m}} - m{v}_{3}^{2}} \right)\frac{{\zeta {{p}_{y}}}}{m} + \frac{{{{{v}}_{3}}}}{m}{{p}_{x}}{{p}_{y}}} \right]} \right. \\ \left. { - \,i\frac{U}{{2{{\mathcal{E}}_{1}}}}\left\{ {\left( {\frac{{{{p}^{2}}}}{{2m}} + m{v}_{3}^{2}} \right)\frac{{{{p}_{x}}}}{m} - \frac{{3\zeta {{{v}}_{3}}}}{{2m}}(p_{x}^{2} - p_{y}^{2})} \right\}} \right), \\ \end{gathered} $$

(33)

$$\begin{gathered} {{D}_{{ty}}}({\mathbf{p}}) = \frac{{e\hbar }}{{2{{\mathcal{E}}_{1}}({\mathbf{p}})\sqrt {\mathcal{E}_{1}^{2}({\mathbf{p}}) - {{U}^{2}}{\text{/}}4} }} \\ \times \left( {\left[ {\left( { - \frac{{{{p}^{2}}}}{{2m}} - m{v}_{3}^{2}} \right)\frac{{\zeta {{p}_{x}}}}{m} + \frac{{{{{v}}_{3}}}}{{2m}}(p_{x}^{2} - p_{y}^{2})} \right]} \right. \\ \left. { - \,i\frac{U}{{2{{\mathcal{E}}_{1}}}}\left\{ {\left( {\frac{{{{p}^{2}}}}{{2m}} + m{v}_{3}^{2}} \right)\frac{{{{p}_{y}}}}{m} - \frac{{3\zeta {{{v}}_{3}}}}{{2m}}{{p}_{x}}{{p}_{y}}} \right\}} \right). \\ \end{gathered} $$

(34)

The components of the velocity operator given by the relation ${{{\mathbf{\hat {v}}}}_{\zeta }}$ = ∂$\hat {H}$/∂${\mathbf{\hat {p}}}$ for the effective 2 × 2 Hamiltonian (3), can be presented by the expressions:

$${{{\hat {v}}}_{{\zeta x}}} = \zeta \left( {\begin{array}{*{20}{c}} 0&{ - \frac{1}{m}(\zeta {{{\hat {p}}}_{x}} - i{{{\hat {p}}}_{y}}) + {{{v}}_{3}}} \\ { - \frac{1}{m}(\zeta {{{\hat {p}}}_{x}} + i{{{\hat {p}}}_{y}}) + {{{v}}_{3}}}&0 \end{array}} \right),$$

(35)

$${{{\hat {v}}}_{{\zeta Y}}} = i\left( {\begin{array}{*{20}{c}} 0&{\frac{1}{m}(\zeta {{{\hat {p}}}_{x}} - i{{{\hat {p}}}_{y}}) + {{{v}}_{3}}} \\ { - \frac{1}{m}(\zeta {{{\hat {p}}}_{x}} + i{{{\hat {p}}}_{y}}) - {{{v}}_{3}}}&0 \end{array}} \right).$$

(36)

The intraband velocity V'(p) at HHG in bilayer AB-stacked graphene is given by the formula:

$$\begin{gathered} {\mathbf{V}}{\kern 1pt} '({\mathbf{p}}) = \left[ {{{{v}}_{3}}{\mathbf{p}} - 3\zeta \frac{{{{{v}}_{3}}p}}{{2m}}{\mathbf{p}}\cos 3\vartheta } \right. \\ \left. { + 3\zeta \frac{{{{{v}}_{3}}{{p}^{3}}}}{{2m}}\sin 3\vartheta \frac{{\partial \vartheta }}{{\partial {\mathbf{p}}}} + 2\frac{{{{{\mathbf{p}}}^{2}}}}{{{{{(2m)}}^{2}}}}} \right]\mathcal{E}_{1}^{{ - 1}}({\mathbf{p}}). \\ \end{gathered} $$

(37)

The author is deeply grateful to prof. H.K. Avetissian for permanent discussions and valuable recommendations. This work was supported by the Science Committee of Ministry of Education, Science, Culture and Sport of RA.