Anomalous Hall effect, magneto-optical properties, and nonlinear optical properties of twisted graphene systems

TL;DR

This paper investigates the anomalous Hall effect, magneto-optical, and nonlinear optical properties of twisted graphene systems, revealing tunable valley polarization effects, giant Faraday/Kerr rotations, and colossal nonlinear responses driven by symmetry breaking and small gaps.

Contribution

It provides a comprehensive analysis of the nonlinear optical responses and magneto-optical effects in twisted graphene systems, highlighting the role of valley polarization and orbital magnetization.

Findings

Valley polarization induces anomalous Hall effect tunable by magnetic fields.

Giant Faraday/Kerr rotations observed in terahertz regime.

Colossal nonlinear optical responses with specific tensor components linked to magnetization.

Abstract

We study the anomalous Hall effect, magneto-optical properties, and nonlinear optical properties of twisted bilayer graphene (TBG) aligned with hexagonal boron nitride (hBN) substrate as well as twisted double bilayer graphene systems. We show that non-vanishing valley polarizations in twisted graphene systems would give rise to anomalous Hall effect which can be tuned by in-plane magnetic fields. The valley polarized states are also associated with giant Faraday/Kerr rotations in the terahertz frequency regime. Moreover, both hBN-aligned TBG and TDBG exhibit colossal nonlinear optical responses by virtue of the inversion-symmetry breaking, the small bandwidth, and the small excitation gaps of the systems. Our calculations indicate that in both systems the nonlinear optical conductivities of the shift currents are on the order of $10^3\,\mu$A/V$^2$; and the second harmonic generation (SHG) susceptibilities are on the order of $10^6\,$pm/V in the terahertz frequency regime. Moreover, in TDBG with $AB\textrm{-}BA$ stacking, we find that a finite orbital magnetization would generate a new component $\sigma^{x}_{xx} $ of the nonlinear photoconductivity tensor; while in $AB$-$AB$ stacked TDBG with vertical electric fields, the valley polarization and orbital magnetization would make significant contributions to the $\sigma^{y}_{xx}$ component of the photoconductivity tensor. These nonlinear photo-conductivities are proportional to the orbital magnetizations of the systems, thus they are expected to exhibit hysteresis behavior in response to out-of-plane magnetic fields.

Figures

Equations

$$H_0^{\mu }=H_{TBG}^{\mu }+H_{mass}$$

$$H_{mass}=\left(\begin{array}{cc}\Delta {\sigma }_z& 0\\ 0& 0\end{array}\right),$$

$$H^{\mu s}=H_0^{\mu }+\mu E_v{\tau }_z+s{E}_s{s}_z,$$

$${\sigma }_{ab}^c={\sigma }_{ab,0}^c+{\sigma }_{ab,z}^c{M}_z.$$

$${\sigma }_{xx,0}^x=-{\sigma }_{yy,0}^x=-{\sigma }_{xy,0}^y=-{\sigma }_{yx,0}^y$$

$${\sigma }_{xx,0}^y={\sigma }_{xy,0}^x={\sigma }_{yx,0}^x=-{\sigma }_{yy,0}^y$$

$${\sigma }_{xx,z}^x=-{\sigma }_{yy,z}^x=-{\sigma }_{xy,z}^y=-{\sigma }_{yx,z}^y$$

$${\sigma }_{xx,z}^y={\sigma }_{xy,z}^x={\sigma }_{yx,z}^x=-{\sigma }_{yy,z}^y,$$

$${\sigma }_{ab}^c(0)=\frac{e^3}{{\omega }^2}\sum _{\Omega =\pm \omega }\sum _{lmn}\int \frac{d\mathbf{k}}{(2\pi {)}^d}\text{Re}[{\varphi }_{ab}(f_l-f_n)\frac{v_{nl}^a{v}_{lm}^b{v}_{mn}^c}{(E_{n\mathbf{k}}-E_{m\mathbf{k}}-i\delta )(E_{n\mathbf{k}}-E_{l\mathbf{k}}+\hslash \Omega -i\delta )}]$$

$${\sigma }_{ab}^c(2\omega )=\frac{e^3}{{\omega }^2}\sum _{\Omega =\pm \omega }\sum _{lmn}\int \frac{d\mathbf{k}}{(2\pi {)}^d}{\varphi }_{ab}(f_l-f_n)\frac{v_{nl}^a{v}_{lm}^b{v}_{mn}^c}{(E_{n\mathbf{k}}-E_{m\mathbf{k}}-2\hslash \Omega -i\delta )(E_{n\mathbf{k}}-E_{l\mathbf{k}}+\hslash \Omega -i\delta )},$$

$$H_{\lambda ,{\lambda }^{\prime }}^{\mu }=\left(\begin{array}{cccc}{h}_1(\mathbf{k})& h_{\lambda }& 0& 0\\ h_{\lambda }^{†}& h_2(\mathbf{k})& U_{\mu }(\mathbf{r})& 0\\ 0& U_{\mu }^{†}(\mathbf{r})& h_3(\mathbf{k})& h_{{\lambda }^{\prime }}\\ 0& 0& h_{{\lambda }^{\prime }}^{†}& h_4(\mathbf{k})\end{array}\right),$$

$$H_{\lambda ,\lambda ;}^{\mu s}=H_{\lambda ,{\lambda }^{\prime }}^{\mu }+\mu E_v{\tau }_z+s{E}_s{s}_z,$$

$${\sigma }_{xx,0}^x(K)=-{\sigma }_{xx,0}^x(K^{\prime })$$

$${\sigma }_{xx,z}^x(K)={\sigma }_{xx,z}^x(K^{\prime })$$

$${\sigma }_{xx,0}^y(K)={\sigma }_{xx,0}^y(K^{\prime })$$

$${\sigma }_{xx,z}^y(K)=-{\sigma }_{xx,z}^y(K^{\prime }),$$

$${\sigma }_{xx,0}^x(K^{\mu })=-{\sigma }_{yy,0}^x(K^{\mu })=-{\sigma }_{xy,0}^y(K^{\mu })=-{\sigma }_{yx,0}^y(K^{\mu })$$

$${\sigma }_{xx,z}^y(K^{\mu })={\sigma }_{xy,z}^x(K^{\mu })={\sigma }_{yx,z}^x(K^{\mu })=-{\sigma }_{yy,z}^y(K^{\mu }),$$

$${\sigma }_{xx}^y(K^{\mu })={\sigma }_{xx,0}^y(K^{\mu })+{\sigma }_{xx,z}^y(K^{\mu })M_z(K^{\mu }),$$

$${\sigma }_{xx}^y(+E_v)=\sum _{\mu =\pm 1}{\sigma }_{xx,0}^y(\mu E_v)+\sum _{\mu =\pm 1}{\sigma }_{xx,z}^y(\mu E_v)M_z^{\mu }$$

$${\sigma }_{xx}^y(-E_v)=\sum _{\mu =\pm 1}{\sigma }_{xx,0}^y(-\mu E_v)-\sum _{\mu =\pm 1}{\sigma }_{xx,z}^y(-\mu E_v)M_z^{-\mu },$$

$$\displaystyle\Big{(}\,\sigma^{y}_{xx}(+E_{v})-\sigma^{y}_{xx}(-E_{v})\,\Big{)}/2\;$$

$$=$$

$$\approx$$

Full text

Twisted bilayer graphene (TBG) has drawn significant attention recently due to the observations of the correlated insulating phases, anomalous Hall effect, and unconventional superconductivity Cao et al. (2018); Sharpe et al. (2019); Choi et al. (2019); Kerelsky et al. (2018); Codecido et al. (2019); Lu et al. (2019); Serlin et al. (2019). At small twist angles, the low-energy states of TBG are characterized by two low-energy bands for each valley and spin degrees of freedom Lopes dos Santos et al. (2007); Bistritzer and MacDonald (2011). Around the “magic angles”, the bandwidths of the low-energy bands become very small, and these nearly flat bands are believed to be responsible for most of the unconventional properties observed in TBG. Numerous theories have been proposed to understand the intriguing phenomena observed in TBG Po et al. (2018a); Yuan and Fu (2018); Koshino et al. (2018); Kang and Vafek (2018); Xu and Balents (2018); Zhang et al. (2019a); Song et al. (2018); Ahn et al. (2018); Po et al. (2018b); Tarnopolsky et al. (2019); Pal et al. (2018); Angeli et al. (2018); Liu et al. (2019a); Lian et al. (2018); Zhang et al. (2019b); Wolf et al. (2019); Sboychakov et al. (2018); Isobe et al. (2018); Xu et al. (2018); Huang et al. (2018); Liu et al. (2018); Rademaker and Mellado (2018); Venderbos and Fernandes (2018); Kang and Vafek (2019); Xie and MacDonald (2018); Jian and Xu (2018); Wu et al. (2018); Bultinck et al. (2019); Kozii et al. (2018); Hazra et al. (2018); Wu (2019); Roy and Jurivcić (2019); Angeli et al. (2019); Wu et al. (2019a); Liu et al. (2019b); Hu et al. (2019); Xie et al. (2019).

Recently unconventional superconducting and correlated insulating behavior, as well as quantum anomalous Hall effect have been observed in twisted double bilayer graphene (TDBG) Shen et al. (2019); Liu et al. (2019c); Cao et al. (2019) and trilayer graphene with hexagonal boron nitride (hBN) substrate Chen et al. (2019a, b), which has stimulated extensive theoretical interests Zhang and Senthil (2019); Lee et al. (2019); Liu et al. (2019d); Koshino (2019); Chebrolu et al. (2019); Haddadi et al. (2019). In particular, it has been proposed that topological flat bands generically exist in twisted double bilayer Lee et al. (2019); Liu et al. (2019d); Koshino (2019) and twisted multilayer graphene systems Liu et al. (2019d), and that the topological flat bands with non-vanishing valley Chern numbers are associated with large and valley-contrasting orbital magnetizatons, which may lead to an orbital ferromagnetic state once the valley symmetry is broken either spontaneously or due to external magnetic fields Liu et al. (2019d). The orbital ferromagnetic states are also believed to exist in TBG aligned with hBN substrate Sharpe et al. (2019); Lu et al. (2019); Xie and MacDonald (2018), in which (quantum) anomalous Hall effect has been observed at 3/4 filling of the flat bands around the magic angle Sharpe et al. (2019); Serlin et al. (2019).

In this context, it is natural to ask how to probe such orbital ferromagnetic states in experiments. Certainly the most salient signature of the orbital magnetic state is the anomalous Hall effect (AHE). In typical magnetic materials, anomalous Hall effect results from the interplay between spin ferromagnetism and spin-orbit coupling (SOC) Nagaosa et al. (2010): time-reversal ($\mathcal{T}$) symmetry is first broken in the spin sector leading to spin magnetizations, then the $\mathcal{T}$ symmetry breaking is transmitted from the spin sector to the orbital sector via SOC. In orbital ferroamgnetic systems, however, $\mathcal{T}$ symmetry is directly broken in the orbital sector, and the AHE does not require any microscopic SOC. One thus expects that the AHE in an orbital ferromagnetic metal would be much more conspicuous than that in a spin ferromagnetic metal, and that the AHE will be quantized as nonzero integers if the orbital magnet is insulating. Such argument also applies to magneto-optical effects. As light is directly coupled with the orbital degrees of freedom of electrons, the magneto-optical effects in an orbital ferromagnet should be much more pronounced than those in a spin ferromagnet. Therefore, we expect that there would be significant magneto-optical responses in both hBN-aligned TBG and TDBG. On the other hand, inversion symmetry is broken in both systems, which allows for nonlinear optical effects Kraut and von Baltz (1979); Sipe and Shkrebtii (2000). It is intriguing to ask whether the orbital magnetization and valley polarization in TBG and TDBG can be probed by nonlinear optical responses. Even without orbital magnetism, the nonlinear optical properties of the twisted graphene systems from structural inversion-symmetry breaking is still an open question, which deserves a comprehensive study.

In this paper, we systematically study the AHE, magneto-optical properties, and nonlinear optical properties of both hBN-aligned TBG and TDBG. We find that in additional to AHE, there are also giant magneto-optical Kerr and Faraday rotations in both systems by virtue of the valley-symmetry breaking and orbital magnetizations. Therefore, we propose that the Faraday and Kerr rotations may be a powerful tool to detect the presence of orbital magnetism in twisted graphene systems. Moreover, we also study the nonlinear optical responses such as the shift current and second harmonic generation (SHG) in hBN-aligned TBG and TDBG. We find that both systems exhibit colossal nonlinear optical responses by virtue of the small bandwidths and the small excitation gaps of the twisted graphene systems. To be specific, our calculations indicate that in either system, the shift-current conductivity $\sigma^{c}_{ab}(0)$ is on the order of $10^{3}\,\mu$A/V2, and the SHG susceptibility $\chi^{c}_{ab}(2\omega)$ is on the order of $10^{6}$ pm/V in the terahertz frequency regime. In TDBG with $AB$-$BA$ stacking, we propose that the non-vanishing valley polarization (orbital magnetization) would generate a new component of the nonlinear photoconductivity $\sigma^{x}_{xx}$; while in $AB$-$AB$ stacked TDBG with vertical electric fields, we find that the orbital magnetization would make significant contributions to the $\sigma^{y}_{xx}$ component. We predict that both of these two components of photoconductivities will exhibit remarkable hysteresis behavior in response to out-of-plane magnetic fields. Therefore, these nonlinear photo-conductivities that are generated by the orbital magnetizations may be considered as strong experimental evidence for the valley polarizations and orbital ferromagnetism in the TDBG systems.

I The TBG system aligned with hexagonal BN substrate

We first study the AHE, magneto-optical properties, and nonlinear optical properties of the hBN-aligned TBG system. We consider the situation that TBG is placed on top of a hBN substrate, and the hBN substrate is aligned with the bottom graphene layer. This is actually the device used in Ref. Sharpe et al., 2019, in which an anomalous Hall conductivity $\sim 2.4e^{2}/h$ has been observed at 3/4 filling of the conduction flat band around the magic angle. The hBN substrate is believed to have two effects on the electronic structures of TBG. First, the alignment of the hBN substrate with the bottom graphene layer would impose a staggered sublattice potential on the bottom layer graphene and break the $C_{2z}$ symmetry, which opens a gap at the Dirac points of the flat bands of the magic-angle TBG. Actually the two flat bands for the $K$ valley acquires nonzero Chern numbers $\pm 1$ ($\mp 1$ for the $K^{\prime}$ valley) once a gap is opened up at the Dirac points. Second, the hBN substrate would generate a new moire pattern, which roughly has the same period as the one generated by the twist of the two graphene layers, but are orthogonal to each other Moon and Koshino (2014). However, the moiré potential generated by the hBN substrate is one order of magnitude weaker than that generated by the twist of the two graphene layers Moon and Koshino (2014); Jung et al. (2014). Therefore, as a leading-order approximation, it is legitimate to neglect the moiré potential generated by the hBN substrate Zhang et al. (2019b); Bultinck et al. (2019). With such an approximation, the effective Hamiltonian for the hBN-aligned TBG system is simplified as

[TABLE]

where $\mu\!=\!\pm 1$ is the valley index, and $H_{TBG}^{\mu}$ represents the continuum Hamiltonian for valley $\mu$ as proposed by Bistrizer and MacDonald Bistritzer and MacDonald (2011), and $H_{mass}$ is the “Dirac mass” term at the bottom layer graphene generated by the hBN substrate, which is expressed as

[TABLE]

where $\Delta$ is the staggered sublattice potential exerted on the bottom graphene layer, which is fixed as $17\,$meV throughout this paper. The details of the continuum Hamiltonian of TBG $H_{TBG}^{\mu}$ is given in Supplementary Material.

I.1 Electronic structures

In order to study the effects of breaking the valley and spin symmetries, we artificially apply valley and spin energy splittings

[TABLE]

where $E_{v}$ and $E_{s}$ are positive real numbers denote the valley and spin splittings, and $\mu\!=\!\pm 1$ and $s\!=\!\pm 1$ represent the valley and spin degrees of freedom respectively. $\tau_{z}$ and $s_{z}$ both denote the third Pauli matrix, and are defined in the valley and spin subspace respectively. The bandstructures at the first magic angle $\theta\!=\!1.05^{\circ}$ are shown in Fig. 1, where the solid blue lines, dashed blue lines, solid red lines, and dashed red lines denote the bandstructures of electrons with {$\mu\!=\!-1$, $s\!=\!-1$}, {$\mu\!=\!-1$, $s\!=\!+1$}, {$\mu\!=\!+1$, $s\!=\!-1$}, and {$\mu\!=\!+1$, $s\!=\!+1$} respectively. The thick dashed gray line denotes the chemical potential for some given valley and spin splittings, which is determined by the charge filling $+3/4$, i.e., filling 7 out of the 8 flat bands including the valley and spin degrees of freedom. When $E_{v}\!=\!0$, $E_{s}\!=\!0$, the bandstructures are spin degenerate for each valley as shown in Fig. 1(a). When $E_{v}\!=\!3\,$meV and $E_{s}\!=\!0$, the bandstructures are shown in Fig. 1(b). Clearly the two valleys have been splitted: the flat bands of the $K$ valley are completely filled, while the conduction band of the $K^{\prime}$ valley is half filled. Later we will show that the AHE, orbital magnetization, and magneto-optical effects would be maximal in such a situation. In Fig. 1(c) we show the bandstructures with $E_{v}\!=\!0$ and $E_{s}\!=\!3\,$meV. We see that the spins have been splitted but the valley symmetry is still preserved. In this situation, both AHE and magneto-optical effects vanish since $\mathcal{T}$ symmetry is broken only in the spin sector, but still preserved in the orbital sector. In Fig. 1(d), we show the bandstructures with $E_{v}\!=\!3\,$meV and $E_{s}\!=\!3\,$meV at +3/4 filling. In this situation, both the spin and valley splittings are strong enough such that both the valley and spin polarizations $\xi_{v}$ and $\xi_{s}$ reach their maximal values with $\xi_{v}=\xi_{s}=1/7$ 111The valley and spin polarizations are defined as: $\xi_{v}=(\sum_{s=\pm 1}\rho_{\mu=-1,s}-\sum_{s=\pm 1}\rho_{\mu=+1,s})/\rho_{t}$, and $\xi_{s}=(\sum_{\mu=\pm 1}\rho_{\mu,s=-1}-\sum_{\mu=\pm 1}\rho_{\mu,s=+1})/\rho_{t}$, where $\rho_{\mu,s}$ is the charge density for valley $\mu$ and spin $s$, and $\rho_{t}=\sum_{\mu,s=\pm 1}\rho_{\mu,s}$ is the total charge density. . Then the system enters a quantum anomalous Hall (QAH) insulating phase with Chern number $-1$. Such a phase has been predicted as the ground state at +3/4 filling of magic-angle TBG based on Hartree-Fock calculations including electrons’ Coulomb interactions Xie and MacDonald (2018); Zhang et al. (2019b).

It is important to note that as a result of the staggered sublattice potential from the hBN substrate, a gap $\sim\!4\,$meV has opened up at the Dirac points $K_{s}$ and $K_{s}^{\prime}$ as clearly shown in Fig. 1. Such a $C_{2z}$ symmetry breaking and the gap opening at the Dirac points are essential in achieving AHE, magneto-optical effect and nonlinear optic effects in the TBG system. If $C_{2z}$ symmetry is preserved, the Hamiltonian for each valley and each spin would have $C_{2z}\mathcal{T}$ symmetry ($\mathcal{T}$ denotes time-reversal), which would enforce the Berry curvature to be zero at every $\mathbf{k}$ point. As a result, both the AHE and magneto-optical effects would be forbidden. The $C_{2z}$ operation also connects the two valleys. If $C_{2z}$ symmetry is preserved, the nonlinear optical response is also prohibited as the contributions from the two valleys would exactly cancel each other.

I.2 Anomalous Hall effect

We take the valley and spin splittings ($E_{v}$ and $E_{s}$ in Eq. (3)) as two free parameters which are varied from 0 to $3\,$meV. Then we study the dependence of AHE, spin/orbital magnetizations, and magneto-optical effects on the valley and spin splittings at +3/4 filling. In Fig. 2(a) we first show the dependence of the anomalous Hall conductivity (in units of $e^{2}/h$) on $E_{v}$ and $E_{s}$. Clearly the anomalous Hall conductivity (AHC) is non-vanishing only if $E_{v}\!>\!0$, otherwise $\mathcal{T}$ symmetry is always preserved in the orbital sector and the contributions from the two valleys always cancel each other. One may notice in Fig. 2(a) that there is a small region in the upper right corner in which $\sigma_{xy}$ is quantized as $-e^{2}/h$. This is the region in which both valley and spin polarizations reach their maximal values $\xi_{v}=\xi_{s}=1/7$, and the system enters a QAH phase with quantized Hall plateau. We note that $E_{v}\!\approx\!E_{s}\!\approx\!2.5\,$meV would be strong enough to approach the QAH phase.

It is also interesting to note that when $E_{v}\!\sim\!3\,$meV and $E_{s}\!\sim\!0$ meV, i.e., in the upper left corner of Fig. 2(a), the AHC is maximal with $\sigma_{xy}\!\sim\!-1.9\,e^{2}/h$, and the system is metallic as shown by the bandstructures in Fig. 2(b). In other words, fixing $E_{v}\!\sim\!3$ meV, $\sigma_{xy}$ would decrease from $-1.9\,e^{2}/h$ to $-e^{2}/h$ as $E_{s}$ increases from 0 to $3\,$meV, and the system would go through a transition from a metal to a QAH insulator. In Ref. Sharpe et al., 2019, the magnitude of the measured AHC $\sim 2.4\,e^{2}/h$ around +3/4 filling, and the system is still metallic. It may suggest that the valley splitting dominates over the spin splitting in their sample, such that the system still stays in a metiallic phase as shown Fig. 1(b), with $|\sigma_{xy}|$ greater than $e^{2}/h$. According to our calculations, increasing the spin splitting may drive the system from the metallic phase to the spin and valley polarized QAH phase for the sample used in Ref. Sharpe et al., 2019. We expect that such a phase transition may be assisted by applying in plane magnetic field which only couples to the spin magnetization.

In Fig. 2(b) we plot the dependence of $\sigma_{xy}$ on the chemical potential with $E_{v}\!=\!3\,$meV and $E_{s}\!=\!0\,$meV. The vertical gray dashed line marks the actual chemical potential at 3/4 filling. We see that as a result of the valley splitting, the bands from the $K$ valley are completely filled, and do not contribute to AHC. However, the conduction band of the $K^{\prime}$ valley is only half filled. By virtue of the giant density of states near the conduction band minimum, the chemical potential is just slightly above the conduction band minimum, thus the AHC contributed by the $K^{\prime}$ valley for each spin species is still close to the quantized value $-e^{2}/h$. This explains the large calculated AHC $\sigma_{xy}\!\sim\!-1.9\,e^{2}/h$ when $E_{v}\!\approx\!3\,$meV and $E_{s}\!\approx\!0$ meV shown in Fig. 2(a).

In Fig. 2(c)-(d) we plot the orbital and spin magnetizations (in units of $\mu_{B}$ per moiré primitive cell) in the parameter space spanned by $E_{v}$ and $E_{s}$. The orbital magnetization can be as large as $\sim-1\,\mu_{B}$ when $E_{v}\!\sim\!3\,$meV and $E_{s}\!\sim\!0\,$meV, and gradually decreases to $\sim 0.1\,\mu_{B}$ as $E_{s}$ increases to $\sim\!3\,$meV. On the other hand, in the QAH phase with $E_{v}\!\approx\!E_{s}\!\sim\!3\,$meV, the spin magnetization is as large as $1\,\mu_{B}$ per moiré cell (see Fig. 2(d)), indicating that the magnetization in the QAH phase would be dominated by the spin component. However, the orbital magnetic order is extremely anisotropic, which breaks a discrete $\mathbb{Z}_{2}$ symmetry; while the spin magnetic order is isotropic due to the absence of atomic SOC, which breaks continuous spin rotational symmetry. Therefore, despite being small in magnitude, the orbital magnetization is expected to be much more robust to thermal fluctuations according to Mermin-Wagner theorem.

I.3 Magneto-optical properties

In Fig. 3(a)-(b) we plot the Faraday and Kerr rotations (denoted by $\theta_{F}$ and $\theta_{K}$) as $E_{v}$ and $E_{s}$ increases from 0 to $3\,$meV (see supp. mat. for details). We consider the case that the incident light is normal to the 2D plane, and we first fix the frequency $\hbar\omega\!=\!0.05\,$eV. As clealry shown in the figures, both $\theta_{F}$ and $\theta_{K}$ vanish when $E_{v}=\!0\!$, and their magnitudes increase with the increase of the valley splittings. When $E_{v}\!\sim\!3$ meV, $\theta_{F}\!\sim\!-0.4^{\circ}$ and $\theta_{K}$ is as large as $9^{\circ}$. In the QAH phase (upper right corner), $\theta_{F}\!\approx\!-0.2^{\circ}$ and $\theta_{K}\!\approx\!5.6^{\circ}$ with the incident light frequency $\hbar\omega\!=\!0.05\,$eV. We note that the calculated Kerr rotation in hBN-aligned TBG is at least an order of magnitude greater than those observed in typical spin ferromagnetic materials Antonov et al. (2004). For example, in $3d$ transition metal compounds such as Fe, Co, Ni, and MnPt3, the maximal Kerr angles are typically on the order of $0.1^{\circ}-1^{\circ}$ over the entire frequency regime Antonov et al. (2004) ; in some magnetic multilayers and heterostructrures such as Co/Pd(Pt) multilayers Antonov et al. (2004), Fe/Au multilayers Antonov et al. (2004), and yttrium-iron garnet thin films Tomita et al. (2006), the measured maximal $\theta_{K}$ is also on the order of $0.1^{\circ}-1^{\circ}$ Antonov et al. (2004). In double-layer CrGeTe3, the Kerr rotation $\theta_{K}\sim 0.0007^{\circ}$ at the light frequency $\hbar\omega\approx 0.35\,$eV Gong et al. (2017); while in single-layer CrI3, the Kerr angle $\theta_{K}\sim 0.3^{\circ}$ at the frequency $\hbar\omega\approx 1.95\,$eV, which has been proposed to arise from excitonic effects Wu et al. (2019b). We see that the calculated Kerr angle in hBN-aligned TBG in the terahertz regime is order of magnitude larger than those of any conventional magnetic materials with small SOC.

In typical spin ferromagnetic materials as mentioned above, the magneto-optical phenomena result from the interplay between spin ferromagnetism and SOC: the SOC transimits the TR symmetry breaking from spin sector to orbital sector, and generates orbital magnetizations. The orbital magnetization and magneto-optical effects would be vanishingly small if SOC amplitude is negligible. However, in most magnetic materials, the effects of SOC are perturbative compared with the bandwidths. Therefore, the observed $\theta_{K}$ are typically very small ($\sim 0.1^{\circ}$) in ferromagnetic transition-metal compounds. On the other hands, in twisted grapehen systems, the Faraday and Kerr rotations directly result from the orbital ferromagnetism as indicated by the significant orbital magnetization shown in Fig. 2(c). The spin and orbital magnetization coexist in the hBN-aligned TBG system and are interwined with each other, but there is no microscopic SOC at the single-particle level. Thus it is expected that the Kerr and Faraday rotations in hBN-aligned TBG (with valley symmetry breaking) would exhibit similar behavior as those in Landau levels Shimano et al. (2013); Morimoto et al. (2009); Ikebe et al. (2010); Volkov and Mikhailov (1985) and in quantum anomalous Hall (QAH) insulators Tse and MacDonald (2011); Okada et al. (2016). However, in hBN-aligned TBG, the difference is that here we do not need external magnetic fields nor any SOC to generate the quantized (anomalous) Hall conductivity; instead, the valley/spin symmetry is expected to be broken spontaneously due to Coulomb interactions.

In Fig. 3(c)-(d) we plot the frequency dependence of $\theta_{F}$ and $\theta_{K}$ for {$E_{v}=3\,$meV, $E_{s}=0$} (blue circles) and {$E_{v}=3\,$meV, $E_{s}=3$ meV} (red diamonds). Both $\theta_{F}$ and $\theta_{K}$ increase dramatically as $\omega$ decreases. In particular, in the QAH phase ($E_{v}\!=\!3\,$meV, $E_{s}\!=\!3$ meV) $\theta_{F}\!=\!0.44^{\circ}$ and $\theta_{K}\!=\!-8.2^{\circ}$ at $\hbar\omega\!=\!0.01\,$eV; while when $E_{v}\!=\!3\,$meV and $E_{s}\!=\!0$, $\theta_{F}\!=\!0.87^{\circ}$, and $\theta_{K}\!=\!-40.04^{\circ}$ for $\hbar\omega\!=\!0.01\,$eV. Actually for an QAH insulator with Chern number $C$, in the limit $\omega\to 0$, $\theta_{F}$ should be quantized as integer multiples of the fine-structure constant $C\alpha\approx C/137\,$rad$\approx C\times 0.42^{\circ}$ Qi et al. (2008); Tse and MacDonald (2011), which is consistent with the results shown in Fig. 3(a). In the QAH phae, the Kerr angle is predicted to be quantized as $\pm\pi/2$ in the low-frequency limit Tse and MacDonald (2011). Indeed our results at lower frequencies ($\hbar\omega<0.01\,$eV) indicate that $\theta_{K}$ tend to approach $90^{\circ}$ as $\omega$ is approaching 0 222In our calculations we have chosen an artificial quasi-particle life time $\tau=\hbar/\delta$, with $\delta=3\,$meV, so that the results are reliable only for $\hbar\omega\gg\delta$. Therefore, the results at frequencies lower than $10\,$meV are not presented as they are no longer reliable when $\hbar\omega\!\sim\!\hbar/\tau$. . It worthwhile to note that both $\theta_{F}$ and $\theta_{K}$ would change signs at $\hbar\omega\sim 0.035\,$eV and $\hbar\omega\sim 0.075\,$eV. This is because the real part of the optical anomalous Hall conductivity $\sigma_{yx}(\omega)$ change sign at $\hbar\omega\approx 0.035\,$eV and $0.075\,$eV, leading to the sign change in the Faraday and Kerr rotations (see Supplemental Material for more details). The sign change of the Faraday and Kerr angles may be an interesting feature which can be easily verified experimentally.

I.4 Nonlinear optical properties

We continue to study the nonlinear optical properties of hBN-aligned TBG. In general, the photocurrent $j^{c}(\omega_{3})$ is related to the time-dependent electric fields of the light via the second-order photoconductivity: $j^{c}(\omega_{3})=\sum_{ab}\sigma^{c}_{ab}(\omega_{3})\,E_{a}(\omega_{1})\,E_{b}(\omega_{2})$, where $a,b,c=x,y$ denotes the spatial directions in Cartesian coordinates 333We consider the case of normal incident light, so that the polarization of light is always within the 2D plane. Therefore, $a,b,c=x,y$, and cannot be along the $z$ direction. , and $\omega_{1}$ and $\omega_{2}$ are the frequencies of the two incident photons, and $\omega_{1}+\omega_{2}=\omega_{3}$ Sipe and Shkrebtii (2000). For monochromatic light, the frequency of the incident photons is fixed as $\pm\omega$, thus there could be two distinct second-order optical processes with $\omega_{3}=0$ or $\omega_{3}=2\omega$, corresponding to the generation of the shift current and the second harmonic generation respectively. Regardless of the microscopic mechanism, the nonlinear optical conductivity tensor $\sigma^{c}_{ab}(0)$ and $\sigma^{c}_{ab}(2\omega)$ have the same properties under symmetry operations, thus they have the same symmetry-allowed components. In particular, one may expand $\sigma^{c}_{ab}$ to the leading order of the orbital magnetization $M_{z}$,

[TABLE]

The only symmetry of hBN-aligned TBG is $C_{3z}$, which restricts $\sigma^{c}_{ab,0}$, $\sigma^{c}_{ab,z}$ and $\sigma^{c}_{ab,zz}$ to the following form

[TABLE]

It turns out that there are only two independent photoconductivities $\sigma^{x}_{xx}$ and $\sigma^{y}_{xx}$ for both shift-current and SHG second-order responses. Each of them include two components: one is independent of the orbital magnetization $M_{z}$, and the other is linear in $M_{z}$. The component that is linear in $M_{z}$ may vary with perpendicular magnetic field, and show hysteresis behavior due to the hysteresis loop of the orbital magnetization. Such a hysteresis behavior will only show up with perpendicular magnetic field since the orbital magnetization is pointing along the $\pm z$ direction.

Microscopically, the shift-current photoconductivity $\sigma^{c}_{ab}(0)$ and the second-harmonic photoconductivity $\sigma^{c}_{ab}(2\omega)$ can be derived using second-order perturbation theory, which are expressed as Zhang et al. (2018)

[TABLE]

and the susceptibility of SHG $\chi^{c}_{ab}(2\omega)\!=\!i\sigma^{c}_{ab}(2\omega)/(2\epsilon_{0}\,\omega)$.

In Fig. 4(a) we plot the frequency dependence of the shift-current photoconductivities at the +3/4 filling of hBN-aligned TBG, where the blue and red markers denote the situations with $E_{v}=E_{s}=0$ and $E_{v}=E_{s}=3\,$meV respectively, and the circles and diamonds represent $\sigma^{x}_{xx}(0)$ and $\sigma^{y}_{xx}(0)$ respectively. We note that the photoconductivities are as large as $\sim\pm 4000\,\mu$A/V2 at $\hbar\omega\!\lessapprox\!0.07$ eV, which is unprecedentedly large. As the frequency increases, the photocondutvities can change sign and decrease to $\sim 10^{2}\,\mu$A/V2 at relatively high frequencies $\hbar\omega\!\sim\!0.2-0.3\,$eV, which is comparable to the calculated shift-current conductivity of bilayer CrI3 in the visible-light frequency regime Zhang et al. (2019c). In Fig. 4(b) we show the imaginary part of the SHG susceptibilities, where the blue and red markers denote situations with $E_{v}=E_{s}=0$ and $E_{v}=E_{s}=3\,$meV (the QAH phase), and the circles and diamonds represent $\chi^{x}_{xx}$ and $\chi^{y}_{xx}$ respectively. At low frequencies $\hbar\omega\!\lessapprox\!0.1\,$ the SHG susceptibilities are extremely large, on the order of $10^{6}\,$pm/V, and they gradually decrease to $\sim\!10^{3}-10^{4}\,$pm/V at higher frequencies. Such colossal SHG susceptibilities are orders of magnitudes larger than those observed in other 2D materials with broken inversion symmetry such as monolayer MoS2 Kumar et al. (2013); Malard et al. (2013), monolayer WSe2 Seyler et al. (2015), WS2 Janisch et al. (2014), and antiferromagnetic bilayer CrI3 Sun et al. (2019). In these 2D materials, the reported SHG susceptibilities are typically on the order of $10^{3}-10^{5}$ pm/V in the visible-light frequency regime. The shift-current and SHG responses at other fillings are on the same order of magnitude as those of 3/4 filling, which we refer the readers to Supp. Mat. for details.

The colossal shift-current and SHG responses shown in Fig. 4 can be interpreted as follows. First, it is straightforward to see from Eq. (6) that the SHG and shift-current responses would be significantly enhanced when the one-photon or two-photon energy ($\omega$ or $2\omega$) is in resonance with some excited electronic states at some $\mathbf{k}$ points. Such a resonant enhancement would be much more pronounced if the energy bands are flat, as the flatness of the bands implies that the all electronic states at different $\mathbf{k}$ points would have the same resonant frequency, thus the resonance effect at different $\mathbf{k}$ points would be summed up. One could imagine having some sets of perfectly flat bands, such as Landau levels, but with broken inversion symmetry. Then if the frequency is in resonance with the Landau-level spacing, one would get enormous shift-current and SHG responses. Such giant nonlinear optical effects have never been observed in Landau levels because there is always inversion symmetry in Landau levels of free 2D electrons’ gas or free Dirac fermions in graphene. In hBN-aligned TBG, however, the low-energy states can be interpreted as pseudo Landau levels Liu et al. (2019a) with broken inversion symmetry due to the presence of the hBN substrate. These bands are roughly flat around the magic angle as shown in Fig. 1, which is expected to contribute to giant nonlinear optical responses once the incident photon frequency is somewhere in resonance with the electronic excitations. The electronic excitations in hBN-aligned TBG turn out to have small gaps, i.e., $\sim\!1\,$meV in the QAH phase at 3/4 filling, and $\sim 15\,$meV at the full filling (see Fig. 1), which indicates that the resonant frequencies in hBN-aligned TBG is very small. The small resonant frequencies would further amplify the giant nonlinear optical responses due to the $1/\omega^{2}$ dependence of the nonlinear photoconductivities (see Eq. (6)).

II Twisted double bilayer graphene

II.1 Model Hamiltonian for TDBG

The continuum Hamiltonian proposed for TBG Bistritzer and MacDonald (2011) can be further generalized to TDBG Chebrolu et al. (2019); Lee et al. (2019); Koshino (2019); Liu et al. (2019d):

[TABLE]

where $h_{l}(\mathbf{k})$ ($l=1,2$) denotes the low-energy effective Hamiltonian for monolayer graphene, i.e., $h_{l}(\mathbf{k})=-v_{F}(\mathbf{k}-\mathbf{K}_{l}^{\mu})\cdot[\mu\sigma_{x},\sigma_{y}]+(l-1)U_{d}/3$, where $\mathbf{K}_{l}^{\mu}$ are the $K$ or $K^{\prime}$ point of the $l$th layer, $v_{F}$ is the bulk Fermi velocity, $U_{d}$ is the vertical electrostatic potential across the double bilayer, and $\mu=\pm 1$ is the valley index. $h_{\lambda}$ is the interlayer coupling of the $AB$ ($\lambda\!=\!+1$) or $BA$ ($\lambda\!=\!-1$) stacked bilayer graphene, and $U_{\mu}(\mathbf{r})$ is the moiré potential for valley $\mu$, which is generated by the twist of the double bilayers. We refer the readers to Supp. Mat. for the explicit expressions of the interlayer coupling $h_{\lambda}$ and the moiré potential $U_{\mu}(\mathbf{r})$. The effective Hamiltonian for TDBG with $AB$-$AB$ and $AB$-$BA$ stackings would correspond to $H_{+1,+1}^{\mu}$ and $H_{+1,-1}^{\mu}$ respectively. Eq. (7) is the Hamiltonian for each valley and spin species. In order to study the effects of breaking the valley and spin symmetries, we apply artificial valley and spin splittings $H_{\lambda,\lambda^{\prime}}^{\mu}$:

[TABLE]

where $E_{v}$ and $E_{s}$ are positive real numbers denote the valley and spin splittings, and $\mu\!=\!\pm 1$ and $s\!=\!\pm 1$ represent the valley and spin degrees of freedom respectively. $\tau_{z}$ and $s_{z}$ both denote the third Pauli matrix, but defined in the valley and spin space respectively. In what follows we will study the dependence of AHC and magneto-optical effects on $E_{v}$ and $E_{s}$.

II.2 Anomalous Hall effect and Magneto-optical properties

II.2.1 $AB$-$BA$ stacking

We first calculate the AHC at zero filling (filling 4 out of the 8 low-energy bands) and +1/2 filling (filling 6 out of the 8 low-energy bands) of $AB$-$BA$ stacked TDBG at the twist angle $\theta\!=\!1.05^{\circ}$. In Fig. 5(a) we show the the AHC of $AB$-$BA$ stacked TDBG at zero filling in the parameter space of $E_{v}$ and $E_{s}$. We see that the AHC is as large as $\sim-3.5\,e^{2}/h$ when the valley splitting $\sim\!3\,meV$, and gradually deceases with the increase of spin splitting. With the Hamiltonian Eq. (8), the TDBG system is always metallic for $0\!\leq\!E_{v}\!\leq\!3\,$meV and $0\!\leq\!E_{s}\!\leq\!3\,$meV, thus the AHC is not quantized. In Fig. 5(b) we show the AHC of $AB$-$BA$ stacked TDBG at +1/2 filling. Again, our calculations show that a small valley splitting $\sim\!3\,$meV would generate a substantial AHE with the AHC $\sim-3.4e^{2}/h$. Certainly the specific value of the AHC is sensitive to the details of the system. However, our calculations indicate that small valley splittings in $AB$-$BA$ stacked TDBG would lead to giant AHE, and such a feature of TDBG should be qualitatively correct. This is because such a property is a direct consequence of the valley Chern numbers and orbital ferroamgnetism of the low-energy bands in TDBG Liu et al. (2019d), and it cannot occur in conventional spin ferromagnetic metals and insulators in which the AHE is generated through SOC as a perturbative effect. In the latter the AHC is proportional to the strength of SOC, which is typically much smaller than the bandwidth. As a result, the AHC is typically much smaller in magnitude than the diagonal conductivities.

The above argument also applies to magneto-optical phenomena. In Fig. 6(a) and (b) we plot the Faraday rotations of $AB$-$BA$ stacked TDBG at $\theta\!=\!1.05^{\circ}$ at zero filling and $+1/2$ filling respectively, with the incident light frequency $\hbar\omega\!=\!0.05\,$eV. For zero filling (Fig. 6(a)), the Faraday angle $\theta_{F}$ is largest around the upper right corner when $E_{s}\!\approx\!E_{v}\!\sim 2-3\,$meV, with the maximal $\theta_{F}\!\sim\!0.5^{\circ}$. For +1/2 filling (Fig. 6(b)), $\theta_{F}$ is largest when $E_{v}\!\approx\!3\,$meV, and $E_{s}\!\sim\!0$ meV, and the maximal $\theta_{F}\!\sim\!0.2^{\circ}$. By virtue of the orbital magnetism and nontrivial valley Chern numbers, small valley splittings $\sim\!2-3\,$meV would be strong enough to generate giant Faraday rotations. In Fig. 6(c) and (d) we plot the Kerr rotations of $AB$-$BA$ stacked TDBG for zero filling and $+1/2$ filling respectively. We find that the valley splittings and orbital magnetization in $AB$-$BA$ stacked TDBG would generate giant Kerr rotations $\sim-10^{\circ}$ at zero filling and $\sim-5^{\circ}$ at +1/2 filling, with the incident light frequency $\hbar\omega=0.05\,$eV.

II.2.2 $AB$-$AB$ stacking

On the other hand, in $AB$-$AB$ stacked TDBG, the orbital magnetization for each valley vanishes as a result of $C_{2x}$ symmetry Liu et al. (2019d), which implies that the AHC and the magneto-optical Kerr/Faraday rotations would vanish as well, as both effects are induced by the orbital magnetization. However, the $C_{2x}$ symmetry for each valley would be broken due to the presence of vertical electric fields, which gives rise to isolated topological flat bands with tuable valley Chern numbers Lee et al. (2019); Koshino (2019); Liu et al. (2019d). The non-zero valley Chern numbers of the topological flat bands are associated with valley contrasting orbital magnetizations, which may lead to substantial AHE and Kerr/Faraday rotations if the valley symmetry is broken.

Ferromagnetic insulating states have been observed in experiments at 1/2 filling of the isolated topological flat bands in $AB$-$AB$ stacked TDBG Cao et al. (2019); Liu et al. (2019c); Shen et al. (2019). Such ferromagnetic insulating states have been proposed to be spin ferromagnetic states Liu et al. (2019c); Lee et al. (2019), which is not expected to exhibit any AHE nor magneto-optical effects due to the negligible SOC in graphene. However, at 1/4 or 3/4 filling of the isolated topological flat bands, it is possible to achieve a QAH state with both valley and spin polarizations. If such a state could be realized, then the system is expected to have significant Faraday and Kerr effects as in the cases of $AB$-$BA$ stacked TDBG and hBN-aligned TBG.

II.3 Nonlinear optical properties

II.3.1 $AB$-$BA$ stacking

We continue to study the nonlinear optical properties of TDBG. Again, before going into the details, we first make symmetry analysis on the nonlinear photo-conductivity tensor. $AB$-$BA$ stacked TDBG has both $C_{2y}$ and $C_{3z}$ symmetries Koshino (2019). The $K$ and $K^{\prime}$ valleys are invariant under $C_{3z}$ operation, but are interchanged with each other under $C_{2y}$ operation. Therefore, for each valley of $AB$-$BA$ stacked TDBG, there is only $C_{3z}$ symmetry, and the symmetry allowed form of the second-order photoconductivity tensor for each valley is already given by Eq. (5). However, in $AB$-$BA$ TDBG, the $C_{2y}$ symmetry would further enforce the photo-conductivities of the $K$ and $K^{\prime}$ valleys to the following form,

[TABLE]

where $\sigma^{x}_{xx,0}$, $\sigma^{x}_{xx,z}$, $\sigma^{y}_{xx,0}$, $\sigma^{y}_{xx,z}$ are defined in Eq. (4). It follows that if the two valleys remain degenerate, $\sigma^{x}_{xx}=\sigma^{x}_{xx,0}+\sigma^{x}_{xx,z}M_{z}$ must vanish, because the contributions from the opposite valleys (with opposite $M_{z}$) would exactly cancel each other. If the valley symmetry is broken, the net orbital magnetization $M_{z}$ would be non-vanishing, leading to nonzero $\sigma^{x}_{xx}\sim\sigma^{x}_{xx,z}M_{z}$. Thus the $\sigma^{x}_{xx}$ component of the nonlinear photoconductivity tensor can be used as a probe to detect the valley symmetry breaking and the associated orbital magnetization in $AB$-$BA$ stacked TDBG. Once the orbital $\mathcal{T}$ symmetry (valley symmetry) is spontaneously broken, the $\sigma^{x}_{xx}$ component would be linearly proportional to the orbital magnetization of the system, which is expected to exhibit hysteresis behavior in response (only) to the out-of-plane magnetic field.

In Fig. (7)(a) we plot the shift-current photoconductivity $\sigma^{x}_{xx}(0)$ at +1/2 filling. The blue circles, red diamonds, and magenta plus signs represent the situations with the valley polarization $\xi_{v}=0$, $\xi_{v}=0.1$, and $\xi_{v}=-0.1$ respectively. Clearly, when the valley symmetry is preserved (valley polarization $\xi_{v}=0$), $\sigma^{x}_{xx}$ identically vanishes at any frequency, as the contributions from the two valleys exactly cancel each other. On the other hand, with $\pm 10\%$ of valley polarizations ($\xi_{v}\!=\!\pm 0.1$), $\sigma^{x}_{xx}$ can be as large as $\pm 10^{3}\,\mu$AV*-2* at relatively low frequencies $\hbar\omega\lessapprox 0.1\,$eV as shown in Fig. (7)(a). Moreover, $\sigma^{x}_{xx}$ are opposite for opposite valley polarizations, because $\sigma^{x}_{xx}$ is proportional to the total orbital magnetization as argued above. At higher frequencies $\hbar\omega\gtrapprox 0.1\,$eV, $\sigma^{x}_{xx}$ gradually decreases to $\sim 10^{2}\mu\,$A V*-2*. In Fig. (7)(b) we show the shift-current photoconductivity $\sigma^{y}_{xx}(\omega=0)$ at 1/2 fillings with the valley polarizations $\xi_{v}=0$, $+0.1$ and $-0.1$. Clearly the valley polarization does not significantly change the value of $\sigma^{y}_{xx}$, which is always on the order of $10^{3}\,\mu$A V*-2* for $\hbar\omega\lessapprox 0.1\,$eV, and decrease to $\sim 10^{2}\,\mu$A V*-2* at higher frequencies. It is worthwhile to note that $\sigma^{y}_{xx}$ are nearly identical for $\xi_{v}\!=\!\pm 0.1$, which is expected according to Eq. (9).

In Fig. 7 (c)-(d) we plot the imaginary part of the SHG susceptibilities $\chi^{x}_{xx}(2\omega)$ and $\chi^{y}_{xx}(2\omega)$ at 1/2 filling of $AB$-$BA$ stacked TDBG. Again, the blue circles, red diamonds, and magenta plus signs denote the cases with 0%, +10% and -10% valley polarizations respectively. We see that $\textrm{Im}\chi^{x}_{xx}(2\omega)$ is giant ($\sim\pm 10^{6}\,$pm/V) when $\xi_{v}=\pm 0.1$ at low frequencies, but vanishes when $\xi_{v}=0$. Again, since $\chi^{x}_{xx}(2\omega)$ is directly proportional to the total orbital magnetization of the system, it has opposite signs for opposite valley polarizations. To the contrary, $\textrm{Im}\chi^{y}_{xx}(2\omega)$ seems to be not sensitive to the valley polarizations, which is always on the order of $10^{6}\,$pm/V for $\hbar\omega\lessapprox\,0.1$eV, and gradually decreases to $10^{3}-10^{4}\,$pm/V at higher frequencies.

II.3.2 $AB$-$AB$ stacking

In $AB$-$AB$ stacked TDBG, there are $C_{3z}$ and $C_{2x}$ symmetries for each valley. As a result, the only non-vanishing photoconductivity components for each valley are $\sigma^{x}_{xx,0}$ and $\sigma^{y}_{xx,z}$ as defined in Eq. (4), i.e.,

[TABLE]

where $\mu=\pm 1$ refers to the valley index, with $K^{-}=K$, and $K^{+}=K^{\prime}$. However, the $C_{2x}$ symmetry further enforces that $M_{z}\!=\!0$ for each valley, which implies $\sigma^{y}_{xx}=\sigma^{y}_{xx,z}\,M_{z}$ must vanish. The $C_{2x}$ symmetry can be broken by vertical electric fields, which would allow for non-vanishing but valley-contrasting orbital magnetizations. If the valley symmetry is further broken either spontaneously by Coulomb interactions or by external magnetic fields, then the orbital magnetization would contribute to $\sigma^{y}_{xx}$ in such a way that $\sigma^{y}_{xx}$ would exhibit hysteresis behavior under out-of-plane magnetic fields. To be specific, in the presence of vertical electric fields, for the valley $K^{\mu}$, $\sigma^{y}_{xx}(K^{\mu})$ is expressed as

[TABLE]

Note that the coefficients $\sigma^{y}_{xx,0}(K^{\mu})$ and $\sigma^{y}_{xx,z}(K^{\mu})$ do not change signs under time-reversal operation, and they are dependent on the valley indices only if the bandstructures and/or chemical potentials of the two valleys are different. Then we consider the situation that the system has nonzero valley splittings $\pm E_{v}$ (see Eq. (8). When the valley splitting is $+E_{v}$, the orbital magnetization of the $K$ and $K^{\prime}$ valleys are denoted as $M_{z}^{-}$ and $M_{z}^{+}$ respectively; if the valley splitting is reversed to $-E_{v}$, then the orbital magnetizations of the $K$ and $K^{\prime}$ valleys would become $-M_{z}^{+}$ and $-M_{z}^{-}$. Plugging these relationships into Eq. (11), one obtains

[TABLE]

where $\sigma^{y}_{xx}(\pm E_{v})$ stands for the total $\sigma^{y}_{xx}$ with $\pm$ valley splittings, summing over the contributions from the two valleys. Subtracting $\sigma^{y}_{xx}(+Ev)$ by $\sigma^{y}_{xx}(-E_{v})$ would eliminate the $\sigma^{y}_{xx,0}$ term, leaving the term that is proportional to the total orbital magnetization of the system, i.e.,

[TABLE]

where $M_{z}^{+}+M_{z}^{-}$ is the total orbital magnetization of the system with $+E_{v}$ valley splitting. In the last line of Eq. (13) we have assumed $\sigma^{y}_{xx}(-E_{v})\approx\sigma^{y}_{xx}(0)\approx\sigma^{y}_{xx}(E_{v})$ for small valley splittings. Similar argument also applies to $\sigma^{x}_{xx}$. When $C_{2x}$ symmetry is broken by vertical electric fields, the $\sigma^{x}_{xx,z}$ parameter would be non-vanishing for each valley, such that the orbital magnetization would also contribute to $\sigma^{x}_{xx}$. One can also extract the term proportional to $M_{z}$ by subtracting $\sigma^{x}_{xx}(-E_{v})$ from $\sigma^{x}_{xx}(+E_{v})$.

In Fig. 8 we show the shift-current response of $AB$-$AB$ stacked TDBG with vertical electrostatic potential drop $U_{d}=0.045\,$eV across the four layers. In Fig. 8(a) we present the shift-current photoconductivity $\sigma^{y}_{xx}(0)$ at 3/4 filling of the isolated conduction flat band at $\theta=1.05^{\circ}$. The blue, red and magenta circles represent the cases with valley splitting $E_{v}\!=\!0$, $+1\,$meV and $-1\,$meV respectively. The black diamonds denote $\sigma^{y}_{xx}(+E_{v})/2-\sigma^{y}_{xx}(-E_{v})/2$ with $E_{v}=1\,$meV, which extracts the orbital-magnetization contribution to $\sigma^{y}_{xx}$ as explained in Eqs. (11)-(13). We see that $\sigma^{y}_{xx}$ is actually dominated by the $\sigma^{y}_{xx,z}M_{z}$ term, which is expected to show remarkable hysteresis loops under out-of-plane magnetic fields. In Fig. 8(b) we plot the dependence of $\sigma^{x}_{xx}(0)$ on the light frequency $\omega$ at the same filling and the same twist angle. Clearly $\sigma^{x}_{xx}(0)$ is not changed too much by the valley splitting $E_{v}$, indicating that the orbital-magnetization contribution ($\sigma^{x}_{xx,z}M_{z}$) plays a minor role in $\sigma^{x}_{xx}(0)$.

III Conclusion

To summarize, we have systematically studied the anomalous Hall effect, magneto-optical properties, and nonlinear optical properties of hBN-aligned TBG and TDBG. We have studied the dependence of AHC on the valley and spin splittings in hBN-aligned TBG, and found that the AHE can be engineered using in-plane magnetic fields. In additional to AHE, we also show that there exists giant magneto-optical effect by virtue of the valley-symmetry breaking and orbital magnetizations. We propose that the Faraday and Kerr rotations may be a powerful tool to detect the presence of orbital magnetism in twisted graphene systems. Moreover, we have also studied the nonlinear optical responses, i.e., the shift current and second harmonic generation in both hBN-aligned TBG and TDBG. Our calculations indicate that both systems exhibit colossal nonlinear optical responses. To be specific, the shift-current photoconductivity $\sigma^{c}_{ab}(0)$ is on the order of $10^{3}\,\mu$A/V2, and the SHG susceptibility $\chi^{c}_{ab}(2\omega)$ is on the order of $10^{6}$ pm/V in the terahertz frequency regime. Such gigantic nonlinear optic responses are by virtue of the inversion symmetry breaking, the presence of the low-energy flat bands, and the small excitation gaps in the twisted graphene systems. In TDBG with $AB$-$BA$ stacking, we propose that the non-vanishing valley polarization and orbital magnetization ($M_{z}$) are associated with $C_{2y}$ crystalline symmetry breaking, which would generate a new component of the nonlinear photoconductivity $\sigma^{x}_{xx}\sim M_{z}$; while in $AB$-$AB$ stacked TDBG with vertical electric fields, the valley polarization and orbital magnetization would make significant contributions to $\sigma^{y}_{xx}$. These new components of photoconductivities generated by the orbital magnetizations would exhibit notable hysteresis behavior in response to out-of-plane magnetic fields, and may be considered as strong and robust experimental evidence for the valley polarized state and orbital magnetism in the TDBG system. Our work is a significant step forward in understanding the optical properties of the twisted graphene systems, and may provide useful guidelines for future experimental works.

Acknowledgements.

J.L. and X.D. acknowledge financial support from the Hong Kong Research Grants Council (Project No. GRF16300918). We thank Shiwei Wu, Hongming Weng, Yang Zhang, Hiroaki Ishizuka, and Liuyan Zhao for invaluable discussions.