High-Order Nonlinear Optical Response of a Twisted Bilayer Graphene

Focusing on the twist angle for the minimal commensurate structure, we perform nonperturbative calculations of electron dynamics in the twisted bilayer graphene (TBG) under intense laser fields. We show that the TBG exhibits enriched high-harmonic generation that cannot occur in monolayer or conventional bilayers. We elucidate the mechanism of these nonlinear responses by analyzing dynamical symmetries, momentum-resolved dynamics, and roles of interlayer coupling. Our results imply nonlinear"Opto-twistronics,"or controlling optical properties of layered materials by artificial twists.

Very recently, the twisted bilayer graphene (TBG) has opened a new avenue in physics of Dirac electrons in condensed matter [51,52].The TBG consists of two sheets of graphene vertically stacked with an artificial twist angle, which enables us to manipulate electronic properties of layered materials [53] as sometimes called "twistronics" [54].The twist angle brings about physical phenomena not present in a monolayer graphene such as superconductivity [51,[55][56][57], Mott-like insulating states [51], to name a few.Microscopic theories [58][59][60][61][62] of the TBG have developed, and many active studies are going on to discover and understand novel phenomena [63].
However, the nonlinear optical response of the TBG, or nonlinear "Opto-twistronics", has not yet been explored well.One theoretical challenge is that numerous electronic bands are involved in the TBG due to the large unit cell of the moiré structure.Recently, the Floquet band engineering has been proposed based on the tightbinding model [65] and the low-energy effective Hamiltonian involving a few bands [66][67][68].Another approach is the perturbation theory for the optical field.In this approach, the circular photogalvanic effect [69], one of the lowest-order nonlinear effects, has been found, but analyzing higher-order effects would become more challenging.
In this Letter, by restricting ourselves to a twist angle resulting in the minimal number of bands, we show that the TBG exhibits higher-order nonlinear responses that cannot happen in monolayer or conventional AA-or ABstacked bilayers.The restriction enables the nonperturbative calculation of electron dynamics in the full number of bands.We explain the nonlinear responses character-istic to the TBG by the dynamical symmetries of the Hamiltonian, and elucidate the mechanism of those responses by the reciprocal-space-resolved analysis and the decomposition of the electric current into the intralayer and interlayer contributions.Model and setup.-Webegin by defining the lattice structure of the TBG that we study in this work.We consider two graphenes, or honeycomb lattices, on top of each other, i.e., the AA-stacked bilayer.We let r (l) i denote each site, where l (= up or low) labels each layer and i does each site within the layer.Thus r (up) i and r (low) i share their x and y components, but differ in their z components: [r The minimal commensurate TBG is obtained by rotating the upper (lower) layer by an angle −θ/2 (θ/2) with θ = 21.79 • about the z-axis as illustrated in Fig. 1(b).Thus each site of the TBG is located at , where R z (ϕ) represents the 3 × 3 rotation matrix about the z-axis by angle ϕ.Here the commensurability means the presence of the exact discrete translation symmetry, and the unit cell contains 28 sites for θ = 21.79 • .For other twist angles, the TBG has incommensurate structures or commensurate ones with larger unit cells.One exception is the 60 • -twist, which gives the conventional AB-stacked bilayer.However, as we will see below, the nonlinear optical responses for this case are similar to the monolayer or the AA-stacked bilayer.
To describe the quantum states of the electrons on the TBG, we adopt the tight-binding model of Refs.[60,70] where |R (l) i denotes the Wannier state at position R (l) i .The transfer integral t(R and its parametrization is taken from Refs.[60,70].By the Fourier transform in the xy-plane, we obtain the reciprocal-lattice representation: is the two-dimensional wave vector and the pair (µ, l) (µ = 1, 2, . . ., 14 and l = up or low) serves as the internal degree of freedom corresponding to each site in the TBG unit cell.
The band structure of our TBG is obtained from the eigenvalues of the 28 × 28 Hamiltonian matrix h µl,νl (k) and shown in Fig. 1(d) (see also Ref. [60]).Throughout this work, we assume the half-filling and set E F = 0. We remark that the Dirac cone at the K point is approximately doubly degenerate besides the spin degeneracy.This degeneracy comes from the Dirac electrons of the upper and lower layers.The interlayer coupling does not affect these Dirac electrons much but causes band splittings away from the K point.Now we introduce the coupling of the TBG to the laser propagating in the z-direction.Considering that the laser wavelength is larger enough than the interatomic distances, we assume that the laser electric field E(t) = (E x (t), E y (t), 0) is homogeneous.Then the coupling energy is given by where e is the elementary charge.The total Hamiltonian in the Fourier representation is given by Ĥtotal We focus on a pulse laser of angular frequency Ω, where f (t) represents a 5-cycle Gaussian envelope function [71] and E 0 approximately gives the peak electricfield amplitude.We set the angular frequency as Ω = 0.3 eV corresponding to a mid-infrared laser widely used in experiments (see, e.g., Ref. [5,40]).The parameter p distinguishes the polarization: p = 0 means the linear polarization along the x direction and p = 1 the circular polarization.
Our simulation protocol is as follows.At the initial time t = t ini ( 0), we take the ground state in which every energy eigenstate with negative (positive) energy is occupied (unoccupied).Since we neglect interactions between electrons, we numerically solve the time-dependent Schrödinger equation for individual occupied state under Ĥtotal (t).To reduce the computational cost, we ignore the time evolution of occupied states well below the Fermi energy (E < E F − 5 Ω) since their contributions to the electric current are small.To analyze the optical response, we consider the electric current and its expectation value High-harmonic generation.-First,we analyze the spectra for the electric current induced by the linearlypolarized laser.Figures 2(a) and (b) show the spectra of the currents parallel (J x ) and perpendicular (J y ) to the electric field, respectively.We observe several peaks at (2m + 1)Ω for J x and at 2mΩ for J y (m ∈ Z).In experiments, the induced current with these harmonic peaks is observed as the high-harmonic generation from the TBG as illustrated in Fig. 1(a).In panels (e) and (h), the nonvanishing harmonics of Jx and Jy for linearly-polarized fields are plotted respectively.In panels (g) and (h), we plot the nonvanishing harmonics of Jx in the circularly-polarized fields at the odd and even orders respectively.In panels (e-h), the solids lines show the eye guides ∝ E n 0 for each n.
The even-order harmonics are characteristic to the TBG and cannot appear in the monolayer or conventional AA-and AB-stacked bilayers [72] that have inversion centers [73] although the interlayer bias can give rise to the even-order harmonics [74,75].The selection rules that J x (J y ) has odd-only (even-only) harmonics are explained by the so-called dynamical symmetry appearing in the limit of t FWHM → ∞ [76,77].Note that our TBG without the laser field has the symmetry under C 2y , i.e., the πrotation about the y-axis (see Fig. 1(b)).In the presence of the linearly-polarized electric field, this symmetry is no longer true, but C 2y combined with the time-translation t → t + T /2 becomes a symmetry transformation.This dynamical symmetry leads to the selection rules together with the fact that J x (J y ) is odd (even) under the transformation (see Supplemental Material for detail [71]).To analyze the amplitude of the n-th harmonic, we define the following quantity: where J(ω) represents the spectrum of some component of electric current.In Figs.2(e) and (f), we plot the harmonic amplitude for n ≤ 8 against the incident field amplitude E 0 .For E 0 ≤ 1 MV/cm, each harmonic amplitude scales as A HH n ∝ E n 0 in line with the perturbation theory [2].On the other hand, in the strong-field regime E 0 ≥ 1 MV/cm, A HH n slightly saturates and deviates from the E n 0 -scaling.In this regime, the harmonic peaks are not very sharp as shown in Figs.1(a) and (b) due to lots of excitations occurring between the bands.
Second, we analyze the case of the circular polarization.The current spectra for J x and J y are shown in Figs.2(c) and (d), in which we find a peculiar selection rule: The harmonics at 3mΩ (m ∈ Z) are prohibited.This selection rule derives from another dynamical symmetry consisting of C 3 , the 120 • -rotation about the z-axis, and the timetranslation t → t+T /3.This dynamical symmetry allows the harmonics only at (3m ± 1)Ω and hence prohibits 3mΩ.This symmetry argument also implies that J (ω = 3m ± 1) are circularly polarized [71,77], and thus we obtain similar harmonic peak heights for J x and J y in Figs.2(c) and (d).The harmonic amplitudes and their saturation behavior are shown in Figs.2(g) and (h).
The peculiar selection rule under the circularlypolarized field is characteristic of the TBG and not present in the monolayer or conventional bilayers.The monolayer and the AA-stacked bilayer have the 6-fold rotational symmetry, and thus the harmonics are allowed only for (6m ± 1)Ω [76].The AB-stacked bilayer also allows only harmonics at (6m ± 1)Ω due to the 3-foldrotation and inversion symmetries.These symmetries forbid the harmonics 3mΩ and 2mΩ, respectively, and the allowed harmonics are only (6m ± 1)Ω.The TBG is less symmetric than the monolayer and conventional bilayers, exhibiting enriched nonlinear optical responses with orders n = 6m ± 2.
Reciprocal-space analysis.-Havingfound the harmonic responses characteristic to the TBG, we now investigate their mechanism.To this end, we look into the harmonic amplitude resolved in the reciprocal space by introducing dω Ω J(k; ω), where J(k; ω) represents some component of the Fourier transform of J (k; t).
Figure 3(a) shows the k-resolved second harmonic amplitude |A HH n=2 (k)| obtained for the circularly-polarized field with E 0 = 0.8 MV/cm.The largest amplitude exists in the vicinity of the K and K points and this tendency is commonly seen for the other harmonic orders n.This observation means that large nonlinear currents are carried by the Dirac electrons (see Fig. 1(d)) consistently with the experimental results showing that the Dirac electrons generate harmonics very efficiently [44,45].
Nevertheless, nonDirac electrons play more significant roles in the second harmonic after the sum over the BZ.To show this, we focus on the 6-fold-rotation sum of the k-resolved harmonics and define S HH n (k) We note that the total harmonic amplitude A HH n is obtained as a weighted sum of S HH n (k). Figure 3(b) shows |S HH n=2 (k)| over the k-space, in which we find that the k points near the K point give small contributions.Indeed the individual Dirac electrons carry large nonlinear currents, but these currents cancel each other very strongly.As a result, the nonDirac electrons in the middle of the BZ give more contributions for the second harmonic.The importance of nonDirac electrons are common with other harmonic orders n = 6m ± 2 that are characteristic to the TBG, whereas the Dirac electrons give dominant contributions for the ordinary harmonics n = 6m ± 1.
The band structure in Fig. 1(d) confirms this interpretation.As noted above, the interlayer coupling, emerging as small band splittings, is more effective away from the K point.Given that the interlayer coupling activates the characteristic harmonics n = 6m±2, they are contributed from the k points away from the K point.
Role of interlayer coupling.-Toelucidate other aspects of the interlayer coupling, we decompose the total electric current into two parts, the intralayer and interlayer contributions, as The definitions of these contributions follow from the fact that the current operator Ĵ (k; t) has a 28×28-matrix representation Ĵ (k; t) = µ,l,ν,l j µl,νl (k) |k; µ, l k; ν, l |.We define the operators Ĵintra and Ĵinter as the l = l and l = l parts of Ĵ (k; t), respectively, and J intra (t) and J inter (t) are their expectation values.Figure 4(a) schematically illustrates Ĵintra and Ĵinter , which are the electric currents accompanied by the intralayer and interlayer hoppings of electrons, respectively.
The intralayer component gives the dominant contribution as shown in Fig. 4(b), which shows the result for the circular polarization with E 0 = 0.8 MV/cm.Since the x and y components are essentially equivalent for the circular polarization, we plot only the x component.
For comparison, we also plot the result for the uncoupled bilayers which are defined by removing all the interlayer hopping, i.e., setting t(R Similarly to the monolayer, the uncoupled bilayers only give the harmonics at n = 6m ± 1.For these harmonics, the difference between the TBG and uncoupled bilayers is quite small, meaning that they are carried by the electrons accelerated within each layer. Remarkably, the dominance of the intralayer current holds also for the harmonics n = 6m ± 2 that are caused by the interlayer coupling.Indeed the interlayer coupling is important and, as shown in Fig. 4(c), there occurs significant charge transfer between the layers including some dc (0Ω) component corresponding to the photogalvanic effect [69].However, the in-plane currents accompanied by the interlayer hopping give less contribution to the total current.Rather, the in-plane currents are contributed more by the intralayer electron hopping, and the interlayer coupling assists them by breaking higher symmetry of the uncoupled bilayers and preventing the harmonic currents from canceling out in the BZ.
Discussions and Conclusions.-Wehave conducted the nonperturbative calculations of the laser-induced electric currents in the minimal commensurate TBG, finding higher-order harmonic responses that are not present in monolayer or conventional bilayers.In contrast to the common harmonics, these new harmonics are carried more by nonDirac electrons and caused by the interplay between the intralayer and interlayer electron hoppings.The selection rules of the harmonics are qualitatively distinct and could be tested within the current optics technology.The enriched harmonics in the TBG offer versatile frequency-conversion channels for future applications.
An important future direction toward nonlinear "Opto-twistronics" is to unravel the dependence on the which has been fixed to θ = • in this work.Qualitative results might be different for smaller angles since there occur some emergent symmetries [78,79].Another direction is to go the deep nonperturbative with even stronger fields.In this regime, one should include relaxation due to, e.g., the interband dephasing [80] and impurity scattering [81,82].We leave these open issues for future

S1. ENVELOPE FUNCTION FOR THE PULSE LASER
The explicit form of envelope function is given by Here t FWHM is the full width at half-maximum of the intensity ∼ E(t) 2 rather than field amplitude |E(t)|.define the cycle of the laser pulse by the ratio t FWHM /T with T ≡ 2π/Ω.We use the 5-cycle pulse (t FWHM /T = 5) to obtain all the data presented in the main text.Note that the time-derivative of f (t) gives a minor correction to E(t) = −dA(t)/dt, since f (t) varies slowly in such a multicycle pulse.Thus E 0 in the main text gives the peak electric-field amplitude approximately.

S2. SIMULATION DETAILS
We take 25 × 25 k-points according to the Monkhorst-Pack method in which the point-group symmetry of the hexagonal lattice is respected.The accessible number of k-points is limited by the computational time, but we have confirmed that the qualitative features of the numerical results do not change by varying the number.
For each k-point, we diagonalize the 28 × 28 Hamiltonian matrix h µl,νl (k), obtaining the eigenstates φ a (k) and the corresponding eigenenergy E a (k) for a = 1, 2, . . ., and 28.We let the median of all the eigenenergies {E a (k)} a,k be E M ( E F ).For each k, we pick up the eigenvalues in the range of The evolution is numerically solved for each (k, a) with a ∈ Λ k .We take the initial time t ini = −15T = −3t FWHM where f (t ini ) 0. With the initial condition ψ a (k, t = t ini ) = φ a (k), we numerically integrate the time-dependent Schrödinger equation i∂ ψ a (k, t)/∂t = h(k + eA(t)) ψ a (k, t) up to the final time t fin = 15T = 3t FWHM .We use the Runge-Kutta method with the time step ∆t = (t fin − t ini )/2 12 .We have confirmed that the results show almost no change with smaller time steps.
The time profile of the electric current J (k; t) is calculated as the sum over each initial state J (k; t) = a∈Λ k J a (k; t).Here, J a (k; t) is calculated as the expectation value of the 28 × 28 current matrix in terms of the solution ψ a (k, t).The total current (density) J (t) is obtained as the average of J (k; t) over the 25 × 25 k-points.Since dissipation is neglected in our calculation, small oscillations remain after the pulse irradiation as illustrated in Fig. S1.This is a common technical problem in dissipation-free models and the standard procedure to eliminate these oscillations is to multiply some window function to the calculated time profile of the electric current (see e.g.Ref. [17]).use a similar technique and describe the procedure by taking, e.g., J x (t) the following.First, we discard each 10% of the data in the and at the of J x (t).Thus we have J x (t) (t ini ≤ t < t fin ) with t fin = 0.8t fin = −t ini .Second, we multiply the hanning window to the raw data J x (t), and obtain Jx (t) ≡ J x (t)w(t), which is shown in Fig. S1.This method safely truncates the unphysical persistent oscillations and the Fourier transform of Jx (t) becomes clear without loosing the important features in the middle of the pulse irradiation.All the spectra in the main text are obtained by this procedure.

S4. DYNAMICAL SYMMETRY AND SELECTION RULE
We prove the selection rules for harmonics based on the Floquet theory and dynamical symmetries.These selection rules and symmetries become exact when the incident field is a continuous wave, i.e., t FWHM → ∞.Thus, in the following, we assume the continuous wave and set f (t) = 1.

A. Linear Polarization
We prove that only odd-order (even-order) harmonics are allowed in J x (J y ) when the incident field is linearly polarized along the x direction.The important lattice symmetry is C 2y , which, at the level of Hamiltonian, means the existence of such a 28 × 28 unitary matrix U 2y that where k ≡ (−k x , k y ).Now we consider the Hamiltonian in the presence of the laser field and introduce the new notation for the Hamiltonian matrix h(k; t) ≡ h(k + eA(t)) for clarity below.We note that the Hamiltonian is periodic h(k; t + T ) = h(k; t).Noting that A(t + T /2) = −A(t), we have the following dynamical symmetry Introducing a similar notation for the current matrix by j(k; t), we have The dynamical symmetry (S4) relates the solutions (i.e., the Floquet states) of the time-dependent Schrödinger equation (TDSE) for k and k to each other.Let us focus on the TDSE for k: The Floquet theorem [83] dictates that the independent solutions can be written as where E F a (k) is the so-called quasienergy and u a (k, t) = u a (k, t + T ) is a periodic function for a = 1, 2 . . ., and 28.Similarly, the TDSE for k is given as for each a with ignoring irrelevant phase factors.Thus the Floquet states of k and k are connected to each other by the time shift and unitary transformation.We assume that ψ F a (k, t) and ψ F a (k , t) are equally populated in the dynamics.This is not exactly the case in general, but the population imbalance typically causes little problem (see, e.g., Ref. [23]).Then the total current consists of the contributions from the pairwise Floquet states: We note, e.g., J x a (k; t + T /2) = J x a (k; t − T /2) due to the periodicity.These properties (S13) and (S14) give the selection rules for the harmonics as follows.For the x component, the n-th harmonic amplitude is given by Since we have obtained these selection rules for each pairwise Floquet states, we have similar rules for the total currents as well.Thus we have proved the selection rules for J x and J y , respectively.

B. Circular Polarization
For the circularly-polarized incident field, we have the selection rule that 3mΩ (m ∈ Z) are prohibited.To prove this, we make a parallel argument for the linear polarization with replacing the C 2y symmetry to the C 3 .Correspondingly, the unitary matrix U 2y is replaced by U 3 satisfying h(k; t) = U 3 h(Rk; t + T /3)U † 3 , (S19) where R is the 2 × 2 matrix representation of the 120

FIG. 1 .
FIG. 1.(a) Schematic illustration of high-harmonic generation in twisted bilayer graphene.(b) Top view of the lattice structure of our TBG.The upper and lower layers rotate respectively by the angles −θ/2 and θ/2 with θ = 21.79 • around a common A site.The parallelogram shows the unit cell involving 28 sites.(c) The central solid hexagon shows the first Brillouin zone (BZ) for the superlattice, and the dotted ones the other BZs.The larger hexagons represent the BZs for the upper and lower graphenes.(d) Electronic band structure around the Fermi energy (set to zero) together with (e) the corresponding density of states [60, 64].

FIG. 2 .
FIG. 2. (a-d) Amplitude spectra for in-plane components of electric current plotted for input electric fields E0 = 0.05 (blue), 0.2 (orange), 0.8 (green), and 1.6 MV/cm (red).The polarization of the input electric field is linear (along x) for (a) and (b) and circular for (c) and (d), and the electric-current component is Jx for (a) and (c) and Jy for (b) and (d).(e-h) Amplitudes of n-th harmonic A HH n plotted against the input field amplitude E0.In panels (e) and (h), the nonvanishing harmonics of Jx and Jy for linearly-polarized fields are plotted respectively.In panels (g) and (h), we plot the nonvanishing harmonics of Jx in the circularly-polarized fields at the odd and even orders respectively.In panels (e-h), the solids lines show the eye guides ∝ E n 0 for each n.
FIG. 3. (a) k-resolved harmonic amplitude |A HH n (k)| for n = 2 over the k-space.(b) Absolute value of the 6-foldrotation sum S HH n (k) for n = 2 (see text for definition).In both panels, we use the extended zone scheme, duplicating the data outside the first BZ.

FIG. 4 .
FIG. 4. (a) Schematic illustration of intralayer and interlayer electric currents in the sideview of the TBG.(b-c) Amplitude spectra of the (b) in-plane (Jx) and (b) out-of-plane (Jz) electric current generated by the circular polarization with E0 = 0.8 MV/cm.Each spectrum represents the total J , intralayer Jintra, and interlayer Jinter currents in the TBG and the total current in the uncoupled bilayers.
FIG. S1.Time profile of Jx(t) obtained by the simulation for the circular polarization with E0 = 0.8 MV/cm.(Orange) Modified data Jx(t) with truncation and multiplication of the hanning window.