High-order harmonic generation from twisted bilayer graphene driven by a midinfrared laser field

Theoretical calculations of high-harmonic generation (HHG) from commensurate twisted bilayer graphene (tBLG) under intense laser fields are performed. The nonlinear electron dynamics in tBLG is considered by solving the Liouville–von Neumann equation for a single-particle density matrix, which combines the full energy bands and momentum matrix elements within the framework of tight-binding approximation. We show that the pump intensity determines the relative magnitude of two components of the harmonic spectrum parallel and perpendicular to driving polarization. The important dependence of HHG on twisted angles and crystal orientations is also presented. Especially in the absence of the relaxation process, the harmonic emission for twisted angles around 10 ° exhibits an evident decrease in efficiency per layer compared to monolayer graphene (MLG), which can be interpreted according to Fermi velocity modification. Our calculation also shows that the relative emission efficiency of different harmonic orders between tBLG and MLG contains redundant information on both the dephasing time and an empirical parameter characterizing the decay of the interlayer electron hopping, thus suggesting an all-optical method for the reconstruction of the two parameters. The reconstruction feasibility is successfully demonstrated by a simple optimization algorithm even if considering the possible experimental uncertainty of both driving pulse parameters and high-harmonic signals. Our results show that HHG spectroscopic characteristics in tBLG might serve as a fingerprint to identify the geometric stacking angle and the electronic interaction between adjacent layers, as well as the strong-field laser induced ultrafast dephasing process.

Figure 1

(a) Top view of the commensurate tBLG lattice structure with a twisted angle of $\theta \approx 21.8^{\circ }$ ($m=1,r=1$). ${\mathbf{a}}_1$ and ${\mathbf{a}}_2$ are the primitive vectors of the bottom layer graphene. ${\mathit{\alpha }}_1$ and ${\mathit{\alpha }}_2$ are the basis vectors of the superlattice. (b) Schematic of the BZ corresponding to panel (a). The red-dashed larger hexagon represents the first BZ of the bottom layer with two inequivalent Dirac points $K,K^{\prime }$. The blue-dashed hexagon, which rotates $\theta \approx 21.8^{\circ }$ counterclockwise relative to the red-dashed one, is the first BZ of the top layer with two inequivalent Dirac points $K_{\theta },K_{\theta }^{\prime }$. The seven smaller hexagons in black are the BZ of superlattice tBLG. $M$ and $M$′ are the middle points of the BZ side $K-K_{\theta }$ and $K^{\prime }-K_{\theta }^{\prime }$, respectively.

Figure 2

The parallel (blue solid) and perpendicular (red dashed) polarization components of the high-harmonic spectrum in tBLG with a twisted angle $\theta \approx 21.8^{\circ }$ under the pulse intensity of (a) $I=8.5\times {10}^8\mathrm{W}/\mathrm{c}{\mathrm{m}}^2$ and (b) $I=3\times {10}^{12}\mathrm{W}/\mathrm{c}{\mathrm{m}}^2$. The Gaussian pump pulse has a central wavelength of 4700 nm, a duration of 90 fs, and is linearly polarized along the $x$ axis in Fig. 1.

Figure 3

The total yield of the high-harmonic spectrum from the third to the ninth order as a function of polarization direction $\xi$, calculated for two components parallel (blue solid) and perpendicular (red dashed) to the driving polarization. The inset shows the $\xi$-dependent parallel component plotted on a linear scale. The laser parameters are the same as in Fig. 2.

Figure 4

The fifth-order average harmonic yield per layer as a function of the polarization angle $\xi$, calculated for the conventional AA (red dashed) and AB stacking (blue dotted) bilayer graphene, as well as the SE-even tBLGs with twist angle $\theta \approx 21.8^{\circ }$ (olive dashed-dotted), $\theta \approx 13.2^{\circ }$ (magenta short-dashed), and $\theta \approx 9.4^{\circ }$ (wine short-dotted), respectively. For comparison, the result for the MLG with the same lattice arrangement as the bottom layer of tBLG is also plotted by a black solid line. Laser parameters are the same as in Fig. 2.

Figure 5

The fifth-order harmonic yield in tBLG ($m=1, r=1$) as a function of polarization angle $\xi$, obtained from the contribution of bottom-layer current (red solid), top-layer current (blue dashed), and interference between two layers (olive short-dashed). Laser parameters are the same as in Fig. 2.

Figure 6

The HHG efficiency ratio $R_N\left(\theta \right)$ as a function of the twist angle $\theta$ under three different pump wavelengths: (a) 3500, (b) 4700, and (c) 5900 nm, calculated for the third (blue circles, HH3), fifth (red squares, HH5), seventh (orange stars, HH7), and ninth (purple pluses, HH9) orders. For each case, the pump intensity is fixed at $I=3\times {10}^{12}\mathrm{W}/\mathrm{c}{\mathrm{m}}^2$ and the pulse duration is set as 90 fs. The dephasing process is not considered.

Figure 7

The Fermi velocity ratio ${\overline{v}}_F\left(\theta \right)$ in tBLGs as a function of the twist angle $\theta$, which is generated for the $r=1$ family by varying $m$ from $m=1$ ($\theta \approx 21.8^{\circ }$) to $m=30$ ($\theta \approx 1.{08}^{\circ }$).

Figure 8

Comparison of the harmonic spectrum in MLG calculated by choosing four different size unit cells, including a two-atom cell associated with the original first BZ in MLG (black solid), together with 28-atom, 76-atom, and 148-atom cells associated with the folded BZ whose shapes are the same as ($m=1, r=1$, red dashed), ($m=2, r=1$, blue dotted), and ($m=3, r=1$, olive dashed-dotted) tBLG, respectively. The pump pulse is linearly polarized along the Γ-$K$ direction, with a central wavelength of 5900 nm, a duration of 90 fs, and a peak intensity of $3\times {10}^{12}\mathrm{W}/\mathrm{c}{\mathrm{m}}^2$. The dephasing time is $T_2=10\mathrm{fs}$.

Figure 9

The HHG efficiency ratio $R_N\left(\theta \right)$ in tBLG ($\theta \approx 21.8^{\circ }$) as a function of dephasing time for the HH3 (blue solid), HH5 (red dashed), HH7 (olive dot), and HH9 (black dashed-dotted), obtained from two different relaxation operators: (a) the gauge-variant decoherence term Eq. (23) and (b) the gauge-covariant form Eq. (24). Laser parameters are the same as in Fig. 8.

Figure 10

HHG efficiency ratio $R_N\left(\theta \right)$ as a function of the interlayer decay length $\lambda$, calculated for the HH3 (blue big dots), HH5 (red squares), HH7 (yellow small dots), and HH9 (purple pluses). Laser parameters are the same as in Fig. 7 and the dephasing time $T_2$ is set to 10 fs.

Figure 11

The optimization process for retrieving the dephasing time $T_2$ and interlayer decay parameter $\lambda$. (a) The 24 sets of local optimal solutions ($T_2, \lambda , {\mathrm{\Delta }}_4$) from different initial values, sorted by the ${\mathrm{\Delta }}_4$ in descending order. (b) Projection of panel (a) on the $T_2\text{−}{\mathrm{\Delta }}_4$ plane. (c) Projection of panel (a) on the $\lambda \text{−}{\mathrm{\Delta }}_4$ plane.

Figure 12

The HHG efficiency ratio $R_N\left(\theta \right)$ as a function of (a) the peak intensity and (b) the pulse duration, calculated for HH3 (blue solid), HH5 (red dashed), HH7 (olive dot), and HH9 (black dashed-dotted). The original $T_2^0=10\mathrm{fs}$ and ${\lambda }^0=0.27\AA$ are used.

Figure 13

The optimization process for retrieving the dephasing time $T_2$ and interlayer decay parameter $\lambda$, in the case of introducing to the original peak intensity and pulse duration 20% uncertainty error. (a) The 24 sets of local optimal solutions ($T_2, \lambda , {\mathrm{\Delta }}_2$) from different initial values, sorted by the ${\mathrm{\Delta }}_2$ in descending order. (b) Projection of panel (a) on the $T_2\text{−}{\mathrm{\Delta }}_2$ plane. (c) Projection of panel (a) on the $\lambda \text{−}{\mathrm{\Delta }}_2$ plane.

Figure 14

Energy bands of different bilayer graphene systems. The black dashed line in each panel represents the MLG energy band for comparison.