Collective Excitations of Quantum Anomalous Hall Ferromagnets in Twisted Bilayer Graphene

TL;DR

This paper develops a microscopic theory for collective excitations in quantum anomalous Hall ferromagnets within twisted bilayer graphene, analyzing spin and valley magnon spectra, stability, and thermal effects.

Contribution

It introduces a detailed microscopic framework for understanding collective excitations and stability in QAH ferromagnets in twisted bilayer graphene, including magnon spectra and thermal effects.

Findings

Spin magnon and valley magnon spectra calculated.

QAH ferromagnet stability verified.

Valley wave mode is gapped, favoring valley polarization.

Abstract

We present a microscopic theory for collective excitations of quantum anomalous Hall ferromagnets (QAHF) in twisted bilayer graphene. We calculate the spin magnon and valley magnon spectra by solving Bethe-Salpeter equations, and verify the stability of QAHF. We extract the spin stiffness from the gapless spin wave dispersion, and estimate the energy cost of a skyrmion-antiskyrmion pair, which is found to be comparable in energy with the Hartree-Fock gap. The valley wave mode is gapped, implying that the valley polarized state is more favorable compared to the valley coherent state. Using a nonlinear sigma model, we estimate the valley ordering temperature, which is considerably reduced from the mean-field transition temperature due to thermal excitations of valley waves.

Equations

$$H_0=\sum _{\mathbf{k},\tau ,s}{\epsilon }_{\mathbf{k},\tau }{c}_{\mathbf{k},\tau ,s}^{†}{c}_{\mathbf{k},\tau ,s}$$

$$H_1=\frac{1}{2A}\sum V_{{\mathbf{k}}_1{\mathbf{k}}_2{\mathbf{k}}_3{\mathbf{k}}_4}^{(\tau {\tau }^{\prime })}{c}_{{\mathbf{k}}_1,\tau ,s}^{†}{c}_{{\mathbf{k}}_2,{\tau }^{\prime },s^{\prime }}^{†}{c}_{{\mathbf{k}}_3,{\tau }^{\prime },s^{\prime }}{c}_{{\mathbf{k}}_4,\tau ,s},$$

$$V_{{\mathbf{k}}_1{\mathbf{k}}_2{\mathbf{k}}_3{\mathbf{k}}_4}^{(\tau {\tau }^{\prime })}=\sum _{\mathbf{q}}V(\mathbf{q})O_{{\mathbf{k}}_1{\mathbf{k}}_4}^{(\tau )}(\mathbf{q})O_{{\mathbf{k}}_2{\mathbf{k}}_3}^{({\tau }^{\prime })}(-\mathbf{q}),$$

$$O_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(\tau )}(\mathbf{q})=\int d\mathbf{r}{e}^{i\mathbf{q}\cdot \mathbf{r}}{\Phi }_{\mathbf{k},\tau }^{\ast }(\mathbf{r}){\Phi }_{{\mathbf{k}}^{\prime },\tau }(\mathbf{r}),$$

$$H_{\text{MF}}=$$

$$E_{\mathbf{k},\tau ,s}=$$

$$-\frac{1}{A}\sum _{{\mathbf{k}}^{\prime }}{V}_{\mathbf{k}{\mathbf{k}}^{\prime }\mathbf{k}{\mathbf{k}}^{\prime }}^{(\tau \tau )}{n}_F(E_{{\mathbf{k}}^{\prime },\tau ,s}),$$

$$|\mathbf{Q}{⟩}_S=\sum _{\mathbf{k}}{z}_{\mathbf{k},\mathbf{Q}}{c}_{\mathbf{k}+\mathbf{Q},+,\downarrow }^{†}{c}_{\mathbf{k},+,\uparrow }|\nu =3⟩$$

$${\mathcal{E}}_S(\mathbf{Q})z_{\mathbf{k},\mathbf{Q}}=\sum _{{\mathbf{k}}^{\prime }}{\mathcal{H}}_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(\mathbf{Q})}{z}_{{\mathbf{k}}^{\prime },\mathbf{Q}},$$

$${\mathcal{H}}_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(\mathbf{Q})}=(E_{\mathbf{k}+\mathbf{Q},+,\downarrow }-E_{\mathbf{k},+,\uparrow }){\delta }_{\mathbf{k},{\mathbf{k}}^{\prime }}-\frac{1}{A}{V}_{{\mathbf{k}}^{\prime }(\mathbf{k}+\mathbf{Q})({\mathbf{k}}^{\prime }+\mathbf{Q})\mathbf{k}}^{(++)},$$

$$\mathcal{L}_{S}=-\int d^{2}\boldsymbol{r}\Big{\{}\frac{\hbar n_{0}}{2}\boldsymbol{\mathcal{A}}[\boldsymbol{m}]\cdot\partial_{t}\boldsymbol{m}+\frac{\rho_{s}}{2}(\boldsymbol{\nabla}\boldsymbol{m})^{2}\Big{\}},$$

$$|\mathbf{Q}{⟩}_V=\sum _{\mathbf{k}}{z}_{\mathbf{k},\mathbf{Q}}{c}_{\mathbf{k}+\mathbf{Q},+,\downarrow }^{†}{c}_{\mathbf{k},-,s}|\nu =3⟩,$$

$${\mathcal{E}}_V(\mathbf{Q})z_{\mathbf{k},\mathbf{Q}}=\sum _{{\mathbf{k}}^{\prime }}{\mathcal{W}}_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(\mathbf{Q})}{z}_{{\mathbf{k}}^{\prime },\mathbf{Q}},$$

$${\mathcal{W}}_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(\mathbf{Q})}=(E_{\mathbf{k}+\mathbf{Q},+,\downarrow }-E_{\mathbf{k},-,s}){\delta }_{\mathbf{k},{\mathbf{k}}^{\prime }}-\frac{1}{A}{V}_{{\mathbf{k}}^{\prime }(\mathbf{k}+\mathbf{Q})({\mathbf{k}}^{\prime }+\mathbf{Q})\mathbf{k}}^{(-+)},$$

$$\displaystyle\mathcal{L}_{V}=-\int d^{2}\boldsymbol{r}\Big{\{}$$

$$\displaystyle+\frac{\rho_{\perp}}{2}[(\boldsymbol{\nabla}\pi_{x})^{2}+(\boldsymbol{\nabla}\pi_{y})^{2}]\Big{\}},$$

$$k_{B}T_{DW}=\frac{2}{\ln(1+\sqrt{2})}J\lambda\approx 2.62\Big{(}\frac{\lambda}{a_{M}}\Big{)}^{2}\Delta_{V},$$

$${\mathcal{H}}_{\tau }=\left(\begin{array}{cc}{h}_{\tau b}(\mathbf{k})& T_{\tau }(\mathbf{r})\\ T_{\tau }^{†}(\mathbf{r})& h_{\tau t}(\mathbf{k})\end{array}\right),$$

$$h_{\tau \ell }(\mathbf{k})=e^{-i\tau \ell \frac{\theta }{2}{\sigma }_z}[\hslash v_F(\mathbf{k}-\tau {\mathbf{\kappa }}_{\ell })\cdot (\tau {\sigma }_x,{\sigma }_y)]+{\Delta }_{\ell }{\sigma }_z/2,$$

$$T_{\tau }(\mathbf{r})=T_{\tau }^{(0)}+e^{-i\tau {\mathbf{b}}_{+}\cdot \mathbf{r}}{T}_{\tau }^{(+1)}+e^{-i\tau {\mathbf{b}}_{-}\cdot \mathbf{r}}{T}_{\tau }^{(-1)}$$

Full text

Introduction.— Twisted bilayer graphene (TBG) near the magic angle hosts a plethora of phenomena, e.g., superconductivity Cao et al. (2018a), correlated insulatorsCao et al. (2018b), nematicity Kerelsky et al. (2019); Choi et al. , large linear-in-temperature resistivityCao et al. (a); Polshyn et al. , quantum anomalous Hall effect (QAHE)Sharpe et al. (2019); Serlin et al. , etc. Due to this richness, TBG and related moiré systems are currently under intense experimental Sharpe et al. (2019); Serlin et al. ; Chen et al. ; Cao et al. (2018a, b); Kerelsky et al. (2019); Choi et al. ; Cao et al. (a); Polshyn et al. ; Yankowitz et al. (2019); Codecido et al. ; Lu et al. ; Tomarken et al. (2019); Xie et al. (2019); Jiang et al. (2019); Shen et al. ; Liu et al. ; Cao et al. (b) and theoretical Xu and Balents (2018); Po et al. (2018); Koshino et al. (2018); Kang and Vafek (2018); Liu et al. (2018); Dodaro et al. (2018); Isobe et al. (2018); You and Vishwanath (2019); Tang et al. (2019); Rademaker and Mellado (2018); Guinea and Walet (2018); González and Stauber (2019); Su and Lin (2018); Ramires and Lado (2018); Tarnopolsky et al. (2019); Ahn et al. (2019); Song et al. (2019); Hejazi et al. (2019); Sherkunov and Betouras (2018); Kang and Vafek (2019); Seo et al. (2019); Lin and Nandkishore ; Peltonen et al. (2018); Lian et al. (2019); Choi and Choi (2018); Wu et al. (2018, 2019a); Wu (2019); Wu and Das Sarma (2019, ); Zhang et al. (2019a); Chittari et al. (2019); Wu et al. (2019b); Xie and MacDonald ; Bultinck et al. ; Zhang et al. (2019b); Liu et al. (2019); Lee et al. ; Wu et al. (2019c); Hazra et al. ; Xie et al. ; Julku et al. ; Hu et al. ; Liu and Dai study. For QAHE, which is the focus of this work, moiré bilayers emerge as a new and clean system Serlin et al. ; Chen et al. to realize Chern insulators at elevated temperatures compared with the magnetic topological insulators Chang et al. (2013).

Moiré superlattices in van der Waals bilayers not only generate nearly flat bands, but also often endow the bands with nontrivial topology. In moiré systems with valley contrast Chern numbers, the enhanced electron Coulomb repulsion effect due to band flattening can spontaneously break the valley degeneracy and therefore, lead to valley polarized states with QAHE Zhang et al. (2019a); Chittari et al. (2019); Wu et al. (2019b); Xie and MacDonald ; Bultinck et al. ; Zhang et al. (2019b); Liu et al. (2019); we term such bulk insulating states as quantum anomalous Hall ferromagnets (QAHF), in analogy with the well known quantum Hall ferromagnets (QHF) Moon et al. (1995); Yang et al. (2006). In pristine TBG, $\hat{C}_{2z}$ symmetry (a two-fold rotation around the out-of-plane axis) combined with time-reversal symmetry forbids Berry curvature. However, this $\hat{C}_{2z}$ symmetry can be explicitly broken when TBG is aligned to the hexagonal boron nitride (hBN) substrate, generating a nonzero valley Chern number. It is in this extrinsic TBG aligned with hBN where the anomalous Hall effect (AHE) Sharpe et al. (2019) and later its quantized version (QAHE) Serlin et al. have been observed at the filling factor $\nu=3$. Here we define $\nu$ as $n/n_{0}$, where $n$ is the electron density, and $n_{0}$ the density for one electron per moiré unit cell.

In this paper, we theoretically study the collective excitations in the TBG QAHF, in order to examine the QAHF stability, and to determine the low-energy excitations that control the transport gap and that limit the ferromagnetic transition temperature. The $\nu=3$ QAHF in extrinsic TBG has two distinct collective excitations, i.e., spin magnons and valley magnons, which involve particle-hole transitions with respectively, a single spin flip and a single valley flip. We calculate the energy spectra separately for the two types of magnons by solving their Bethe-Salpeter equations. The calculated excitation spectra indicate that the TBG QAHF is generally robust against small particle-hole fluctuations when the bulk Hartree-Fock gap ($\Delta_{\text{HF}}$) is finite. The spin magnon spectrum has a gapless spin wave mode, which is the Goldstone mode due to the spontaneous breaking of the spin SU(2) symmetry in the $\nu=3$ QAHF. We extract spin stiffness from the long-wavelength spin wave dispersion, and estimate the skyrmion energy. We find that the energy $\Delta_{\text{pair}}$ for a pair of skyrmion and antiskyrmion in the TBG QAHF is comparable in energy with $\Delta_{\text{HF}}$, and either $\Delta_{\text{pair}}$ or $\Delta_{\text{HF}}$ can be the lowest charged excitation gap depending on details of the system.

In a two-dimensional system such as TBG with spin SU(2) symmetry, the spin ordering temperature vanishes according to the Mermin-Wagner theorem. However, QAHE in TBG can arise purely from an orbital effect, e.g., valley polarization. An important distinction between spin and valley is that there is only a valley U(1) symmetry in TBG in contrast to the spin SU(2) symmetry. The $\nu=3$ QAHF preserves the valley U(1) symmetry, but breaks the discrete time-reversal symmetry, which allows a finite valley ordering temperature $T_{V}$. We estimate $T_{V}$ based on the fully gapped valley magnon spectrum, and find that $T_{V}$ is reduced from the mean-field transition temperature due to thermal excitations of valley waves, which provides an explanation for the experimentally observed hierarchy that the transport energy gap of the TBG QAHF is larger than the corresponding Curie temperature Serlin et al. .

Ferromagnetism.— We calculate the moiré band structure of TBG using the continuum Hamiltonian Bistritzer and MacDonald (2011), with details given in the Supplemental Material (SM) SM . We use parameter $\Delta_{b}$ ($\Delta_{t}$) to describe the sublattice potential difference in the bottom (top) graphene layer, and take $(\Delta_{b},\Delta_{t})=(30,0)$ meV Hunt et al. (2013) in order to simulate the experimental situation Sharpe et al. (2019); Serlin et al. where TBG is in close alignment to one of the two (either top or bottom) encapsulating hBN layers. The corresponding moiré band structure in $+K$-valley at twist angle $\theta=1.15^{\circ}$ is shown in Fig. 1, where the first moiré conduction and valence bands are separated by an energy gap about 4 meV (opened up by $\Delta_{b}$), and respectively carry a Chern number $\mathcal{C}$ of $+1$ and $-1$. Because of time reversal symmetry, the first moiré conduction (valence) band in $-K$ valley has a $\mathcal{C}$ value of $-1$($+1$).

We study a minimal interacting model by retaining only the first moiré conduction band states, assuming that all valence band states are filled. The projected Hamiltonian $H$ has the single-particle term $H_{0}$ and the interacting term $H_{1}$,

[TABLE]

where $c^{\dagger}_{\boldsymbol{k},\tau,s}$, $\varepsilon_{\boldsymbol{k},\tau}$ and $\Phi_{\boldsymbol{k},\tau}$ are respectively the fermion creation operation, moiré band energy, and wave function of the first conduction band state with spin label $s$, valley index $\tau$ and momentum $\boldsymbol{k}$. Due to the time reversal symmetry, $\varepsilon_{\boldsymbol{k},\tau}=\varepsilon_{-\boldsymbol{k},-\tau}$ and $\Phi_{\boldsymbol{k},\tau}$=$\Phi^{*}_{-\boldsymbol{k},-\tau}$, where $\boldsymbol{k}$ is measured relative to the moiré Brillouin zone center $\bar{\Gamma}$ point. In $H_{1}$, $A$ is the system area, $O_{\boldsymbol{k}\boldsymbol{k}^{\prime}}^{(\tau)}(\boldsymbol{q})$ is the density matrix element, and $V(\boldsymbol{q})$ is the screened Coulomb potential $2\pi e^{2}\tanh(qd)/(\epsilon q)$, where $\epsilon$ is the effective dielectric constant, and $d$ is the vertical distance between TBG and top(bottom) metallic gates. We take $d$ to be 40 nm as in the experiment of Ref. Serlin et al. , and $\epsilon$ as a free parameter since screening in TBG can be quite complicated. The dielectric screening from the encapsulating hBN should set a lower bound on $\epsilon$, leading to $\epsilon>5$ in our model. $\epsilon$ also effectively controls the ratio between interaction and bandwidth. In TBG the bandwidth near the magic angle is not exactly known experimentally, which is another good reason to take $\epsilon$ as a free parameter.

Hamiltonian $H$ has spin SU(2) and valley U(1) symmetry. We use the Hartree-Fock (HF) approximation and assume that the valley U(1) symmetry is preserved, but allow spin and valley polarization, which leads to the following mean-field Hamiltonian

[TABLE]

where the quasiparticle energy $E_{\boldsymbol{k},\tau,s}$ includes moiré band energy, and Hartree as well as Fock self energies, and $n_{F}$ is the Fermi-Dirac occupation number.

We focus on integer filling factors $\nu=$1, 2 and 3, and make a zero-temperature ($T=0$) ground state ansatz that $\nu$ out of the 4 first moiré conduction bands (including spin and valley degeneracies) are filled, while the remaining $4-\nu$ bands are empty. At $\nu=1$ and 3, the ansatz leads to maximally spin and valley polarized states, which are QAHF and also exact eigenstates of the Hamiltonian $H$. At $\nu=2$, this ansatz generates two distinct types of states, namely, valley polarized state with QAHE and valley unpolarized state without QAHE, which are energetically degenerate at this particular filling, but it is conceivable that a short-range atomic scale interaction (not included in our Hamiltonian $H$) may break this degeneracy. We calculate the $T=0$ HF energy gap $\Delta_{\text{HF}}$ between empty and occupied bands, as shown in Fig. 2(a). A positive $\Delta_{\text{HF}}$ indicates the above ansatz is a good candidate for ground states at least in the HF approximation. As expected, $\Delta_{\text{HF}}$ decreases with increasing dielectric constant $\epsilon$ because of the decreasing interaction strength. $\Delta_{\text{HF}}$ has a strong filling factor dependence, mainly because the Hartree self energy varies strongly with the electron density Guinea and Walet (2018). The gap $\Delta_{\text{HF}}$ at $\nu=3$ is smaller compared to those at $\nu=1$ and 2 for small $\epsilon$, but this order is reversed for large $\epsilon$. By fitting to the experimental $\nu=3$ gap ($\sim$ 2 meV) reported in Ref. Serlin et al. , we estimate $\epsilon$ to be about 30 in our model. With this value of $\epsilon$, we find that the $\nu=1$ and 2 states are not fully gapped in contrast to $\nu=3$, which is consistent with experimental findings in Ref. Serlin et al. . Therefore, our minimal model does capture the essential experimental phenomenology Serlin et al. provided $\epsilon$ is tuned to simulate screening of Coulomb interaction, most likely by all the other moiré bands neglected in our theory.

We show the calculated mean-field ferromagnetic transition temperature $T_{\text{MF}}$ at $\nu=3$ in Fig. 2(b). $T_{\text{MF}}$($\epsilon=30$) is about 22 K, which is larger than the experimental Curie temperature ($\sim$9 K) Serlin et al. . We argue that this discrepancy is due to valley wave excitations, which limit the valley ordering temperature, as will be discussed in the following. The anomalous Hall conductivity $\sigma_{yx}$ at $\nu=3$ is quantized to $e^{2}/h$ within a 0.3% accuracy up to $T=3\text{K}$ as shown in Fig. 2(b). We numerically find that $T_{\text{MF}}$ marks a first-order transition between phases with and without spin-valley polarization, which leads to a jump in $\sigma_{yx}$ at $T_{\text{MF}}$ [Fig. 2(b)]. Remarkably, the experimental anomalous Hall resistance $R_{xy}$ in Ref. Serlin et al., also displays a sizable jump near the Curie temperature.

Spin wave.— We examine the stability of the QAHF by studying the collective excitation spectrum. The spin magnon state at $\nu=3$ can be parametrized as follows

[TABLE]

where $|\nu=3\rangle$ is the QAHF state in which only the valley $+K$ and spin $\downarrow$ band is empty, $z_{\boldsymbol{k},\boldsymbol{Q}}$ are variational parameters, and $\boldsymbol{Q}$ defined within the first moiré Brilouin zone is the momentum of the magnon. In the magnon state $|\boldsymbol{Q}\rangle_{S}$, we make a single spin flip from the occupied spin $\uparrow$ band to unoccupied spin $\downarrow$ band within the same $+K$ valley. Variation of the magnon energy with respect to $z_{\boldsymbol{k},\boldsymbol{Q}}$ leads to the following eigenvalue problem

[TABLE]

where the first part in $\mathcal{H}_{\boldsymbol{k}\boldsymbol{k}^{\prime}}^{(\boldsymbol{Q})}$ is the quasiparticle energy cost of the particle-hole transition, and the second part represents the electron-hole attraction. Equation (4) is typically called the Bethe-Salpeter equation, representing repeated electron-hole interactions (“ladder diagrams”), in the context of excitons in semiconductors; here it gives rise to the spin wave spectrum. We note that $\mathcal{H}_{\boldsymbol{k}\boldsymbol{k}^{\prime}}^{(\boldsymbol{Q})}$ is not gauge invariant (except at $\boldsymbol{Q}$=0) due to the phase ambiguity of the wave function. However, only closed loops in the momentum space appear in the characteristic polynomial of $\mathcal{H}_{\boldsymbol{k}\boldsymbol{k}^{\prime}}^{(\boldsymbol{Q})}$, making its eigenvalues gauge invariant; products of wave function overlap along the closed loops encode information of Berry curvature and quantum geometry Srivastava and Imamoğlu (2015).

We numerically solve Eq. (4), and show the spin excitation spectrum in Fig. 3. The lowest energy mode (spin wave) is gapless at $\boldsymbol{Q}=0$, which is expected from the Goldstone’s theorem, as the continuous spin SU(2) symmetry is spontaneously broken in the QAHF. Because of the spin SU(2) symmetry, the spin lowering operator $\sum_{\boldsymbol{k}}c_{\boldsymbol{k},\tau,\downarrow}^{\dagger}c_{\boldsymbol{k},\tau,\uparrow}$ commutes with the Hamiltonian $H$. Therefore, $z_{\boldsymbol{k},\boldsymbol{Q}}=1$ for any $\boldsymbol{k}$ is an exact zero-energy solution to Eq. (4) at $\boldsymbol{Q}=0$. The overall spin excitation spectrum is nonnegative in the parameter space that we have explored ($\epsilon$ up to 30), showing the stability of the QAHF at $\nu=3$ against spin wave excitations.

The spin wave mode can be phenomenologically described using an $O$(3) nonlinear sigma model Girvin

[TABLE]

where the unit vector $\boldsymbol{m}$ represents the local spin polarization, $\mathcal{A}[\boldsymbol{m}]$ is the effective spin gauge field defined by $\boldsymbol{\nabla}_{\boldsymbol{m}}\times\mathcal{A}[\boldsymbol{m}]=\boldsymbol{m}$, and $\rho_{s}$ the spin stiffness. We estimate $\rho_{s}$ by fitting the numerical spin wave spectrum around $\boldsymbol{Q}=0$ shown in Fig. 3 to the analytical spin wave dispersion $\mathcal{E}_{SW}=(2\rho_{s}/n_{0})\boldsymbol{Q}^{2}$ given by Eq. (5). Besides spin waves, the Lagrangian $\mathcal{L}_{S}$ also supports skyrmion excitations, which are expected to be charged in the case of QAHF, similar to the QHF case Moon et al. (1995). A pair of skyrmion and antiskyrmion has a total energy cost of $\Delta_{\text{pair}}=8\pi\rho_{s}$. We calculate $\Delta_{\text{pair}}$ using $\rho_{s}$ estimated above, and find that $\Delta_{\text{pair}}$ is comparable in magnitude to $\Delta_{\text{HF}}$, but the former is larger at $\nu=3$, as shown in Fig. 2(a). We find that the same ($\Delta_{\text{pair}}>\Delta_{\text{HF}}$) is also true at $\nu=1$ and 2 for spin maximally polarized states. By comparison, $\Delta_{\text{pair}}$ is half of $\Delta_{\text{HF}}$ for the $\nu=1$ quantum Hall ferromagnet in the lowest Landau level with Coulomb interaction Moon et al. (1995). An important difference here with the lowest Landau level is that electron density in the moiré band is spatially nonuniform with modulation within the moiré unit cell, and both the Hartree and Fock self energies modify the moiré bandwidth. Nevertheless, we find that $\Delta_{\text{pair}}$ can be tuned to be smaller than $\Delta_{\text{HF}}$ in TBG by taking both $\Delta_{b}$ and $\Delta_{t}$ to be finite (30 meV), which can be realized when both the top and bottom encapsulating hBN layers are in close alignment to TBG (see SM SM for details). Therefore, we conclude that the lowest charged excitation is determined by either $\Delta_{\text{HF}}$ or $\Delta_{\text{pair}}$, depending on system details.

Valley wave.— Besides spin magnon states, there are also valley magnon states with a single valley flip

[TABLE]

where $s$ can be either $\uparrow$ or $\downarrow$, since both spin components in $-K$ valley are fully occupied in $|\nu=3\rangle$. States $|\boldsymbol{Q}\rangle_{V}$ with $s=\uparrow$ and $\downarrow$ are energetically degenerate for the Hamiltonian $H$, because it actually has an enlarged spin SU(2)$\times$SU(2) symmetry (independent spin rotation within each valley). The corresponding Bethe-Salpeter equation is given by

[TABLE]

which leads to the valley excitation spectrum in Fig. 4. In contrast to the spin excitation spectrum, the lowest valley excitation mode (valley wave) is gapped, consistent with the fact that there is no continuous symmetry broken in the valley pseudospin space. The positive-energy valley wave indicates the robustness of $\nu=3$ QAHF against small variation in the valley space, which implies that the valley polarized state is energetically more favorable than the valley coherent state Bultinck et al. . The valley wave can again be described by a nonlinear sigma model but with an Ising anisotropy

[TABLE]

where the unit vector $\boldsymbol{\pi}$ represents the local valley polarization ($\pi_{z}$ for valley Ising order and $\pi_{x,y}$ for valley coherent order), $u>0$ captures the Ising anisotropy, $\rho_{z,\perp}$ are anisotropic valley stiffness, and other terms are similar to those in Eq. (5). The analytical valley wave dispersion is $\mathcal{E}_{VW}=\Delta_{V}+(2\rho_{\perp}/n_{0})\boldsymbol{Q}^{2}$, where $\Delta_{V}=4u/n_{0}$. Therefore, we can estimate $u$ and $\rho_{\perp}$ using the numerical valley excitation spectrum in Fig. 4.

Because of the Ising anisotropy, there can be valley domain excitations. We make a domain wall ansatz $(\pi_{x},\pi_{y},\pi_{z})=[\text{sech}(x/\lambda),0,\tanh(x/\lambda)]$, and its energy cost is minimized by taking the domain wall width $\lambda$ to be $\sqrt{(\rho_{\perp}+2\rho_{z})/(6u)}$. The domain wall energy per length is then $J=4u\lambda$. We note that this domain wall separates regions with opposite Chern numbers, and binds one-dimensional chiral electronic states. The valley Ising ordering temperature limited by the proliferation of domain walls can be estimated to be Chaikin and Lubensky ; Li et al. (2014)

[TABLE]

where $a_{M}$ is the moiré period. $\Delta_{V}$ can be directly extracted from the valley wave spectrum, but $\lambda$ can not because $\mathcal{E}_{VW}$ has no dependence on $\rho_{z}$. Since $a_{M}$ is the lattice scale in our problem, we argue that the domain wall width $\lambda$ is larger than $a_{M}$, and therefore, we estimate that $k_{B}T_{DW}>2.62\Delta_{V}$. On the other hand, valley waves are already thermally excited when $k_{B}T$ exceeds $\Delta_{V}$. Therefore, we conclude that the valley ordering temperature $T_{V}$ is mostly limited by valley waves instead of domain walls, and estimate $k_{B}T_{V}$ from the valley wave minimum energy. The resulting $T_{V}$ is shown in Fig. 2(b), which is below the mean-field transition temperature $T_{\text{MF}}$. For a zero-temperature charged excitation gap of 2 meV, we find a corresponding $T_{V}$ of about 11 K, which compares well with the experimental Curie temperature Serlin et al. . Although this good quantitative agreement with experiment might be a coincidence, our work establishes the emergent TBG QAHF to be likely a valley Ising ordered state. Regarding the finite jump in the experimentally measured $R_{xy}$ near the transition temperature Serlin et al. , we provide a possible theoretical scenario that the interplay between the continuous spin order parameter $\boldsymbol{m}$ and the Ising degree-of-freedom $\pi_{z}$ through higher-order coupling terms (not included in $\mathcal{L}_{S}$ and $\mathcal{L}_{V}$) could change the finite-temperature phase transition from second order to first order Kato et al. (2010).

Discussion.— In summary, we present a microscopic theory for spin and valley waves of QAHF in TBG, and demonstrate that the excitation spectra provide important information about the stability of mean-field state, the transport energy gap, and the valley ordering temperature. We find that TBG QAHF is robust, provided that other effects such as disorder can be neglected. Besides ferromagnetism, flat moiré bands can host a rich set of broken symmetry states. Our theory can be generalized to study collective excitations of other broken symmetry states in moiré materials.

F. W. thanks A. Young, M. Zaletel, N. Bultinck, S. Chatterjee and I. Martin for discussions. This work was initiated at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611. We acknowledge support by the Laboratory for Physical Sciences.

Note added. While this paper was being written, three related arXiv preprints Repellin et al. ; Alavirad and Sau ; Chatterjee et al. appeared. In this paper we addressed valley ordering temperature limited by valley wave excitations, which has not been studied previously in TBG to our knowledge.

S1 Moiré Hamiltonian

In twisted bilayer graphene (TBG) with a small twist angle $\theta$, the continuum moiré Hamiltonian Bistritzer and MacDonald (2011) is given by

[TABLE]

where $\boldsymbol{r}$ and $\boldsymbol{k}$ are respectively position and momentum operators, and $\tau=\pm$ is the valley index. $h_{\tau b}$ and $h_{\tau t}$ are the Dirac Hamiltonians of the bottom ($b$) and top ($t$) layers:

[TABLE]

where $\ell$ is $+1$ ($-1$) for the $b$ ($t$) layer, $v_{F}$ is the bare Dirac velocity($\sim 10^{6}$ m/s), and $\sigma_{x,y,z}$ are Pauli matrices in the sublattice space. Because of the relative rotation between the two layers, the Dirac cone in layer $\ell$ and valley $\tau$ is shifted to momentum $\tau\boldsymbol{\kappa}_{\ell}=\tau[4\pi/(3a_{M})](-\sqrt{3}/2,-\ell/2)$, which is measured relative to the center of the moiré Brillouin zone [illustrated in Fig. S1(a)]. Here $a_{M}$ is the moiré period given by $a_{0}/[2\sin(\theta/2)]$, where $a_{0}$ is the lattice constant of monolayer graphene. The term $\Delta_{\ell}\sigma_{z}/2$ in Eq. (S2) describes the sublattice potential difference in layer $\ell$. In pristine TBG, both $\Delta_{b}$ and $\Delta_{t}$ vanish. We assume that $\Delta_{b}$ ($\Delta_{t}$) can be induced when TBG is in close alignment to the bottom (top) hBN layers. Under this assumption, $\Delta_{b}$ and $\Delta_{t}$ can be independently controlled.

The interlayer tunneling terms $T_{\tau}(\boldsymbol{r})$ vary in space, following the periodicity of the moiré pattern:

[TABLE]

where $\boldsymbol{b}_{\pm}$ are moiré reciprocal lattice vectors given by $[4\pi/(\sqrt{3}a_{M})](\pm 1/2,\sqrt{3}/2)$ and $T_{\tau}^{(j)}=w_{AA}\sigma_{0}+w_{AB}\cos(2\pi j/3)\sigma_{x}+\tau w_{AB}\sin(2\pi j/3)\sigma_{y}$. Here $w_{AA}$ and $w_{AB}$ are two parameters that respectively determine the tunneling in AA and AB/BA regions of the moiré pattern. We take $w_{AA}=90$ meV and $w_{AB}=117$ meV Wu (2019); Wu et al. (2019a).

The moiré Hamiltonian $\mathcal{H}_{\tau}$ is spin independent, and respects the spinless time-reversal symmetry that relates the two valleys. In the absence of $\Delta_{b}$ and $\Delta_{t}$ , $\mathcal{H}_{\tau}$ builds in the $D_{6}$ point group symmetry that is generated by a sixfold rotation $\hat{C}_{6z}$ around the $\hat{z}$ axis and a twofold rotation $\hat{C}_{2x}$ around the $\hat{x}$ axis. Here the $\hat{C}_{2x}$ operation swaps the two layers. For generic values of $(\Delta_{b},\Delta_{t})$, the point group is reduced to $C_{3}$, since the twofold rotation $\hat{C}_{2z}$, which exchanges the two sublattices within each layer, is broken by the sublattice-dependent potentials $\Delta_{b,t}$. In the special case when $\Delta_{b}=\Delta_{t}$ ($\Delta_{b}=-\Delta_{t}$), the point group is $D_{3}$ that is generated by $\hat{C}_{3z}$ and $\hat{C}_{2y}$ ($\hat{C}_{2x}$).

The Dirac points located at the moiré Brillouin zone corners $\bar{K}$ and $\bar{K}^{\prime}$ are gaped out when the $\hat{C}_{2z}$ symmetry is broken by $\Delta_{b,t}$. We denote the Dirac gap at $\bar{K}$ ($\bar{K}^{\prime}$) in $+K$ valley as $\Delta_{\bar{K}}$ ($\Delta_{\bar{K}^{\prime}}$). If the interlayer tunnelings in TBG were absent, we would find that $\Delta_{\bar{K}}=\Delta_{b}$ and $\Delta_{\bar{K}^{\prime}}=\Delta_{t}$. Due to the interlayer tunnelings, wave functions of TBG are different from those of two decoupled monolayer graphene, and therefore, $\Delta_{\bar{K}}$ ($\Delta_{\bar{K}^{\prime}}$) deviates from $\Delta_{b}$ ($\Delta_{t}$). In Fig. S1(c), we plot $\Delta_{\bar{K}}$ and $\Delta_{\bar{K}^{\prime}}$ as a function of $\Delta_{b}$, while $\Delta_{t}$ is set to 0. In this case $\Delta_{\bar{K}}$ and $\Delta_{\bar{K}^{\prime}}$ take different values, and both are finite, but smaller compared to $\Delta_{b}$. Fig. S1(d) is a similar plot but with $\Delta_{b}=\Delta_{t}$, and therefore, $\Delta_{\bar{K}}=\Delta_{\bar{K}^{\prime}}$.

We obtain the wave function $\Phi_{\boldsymbol{k},\tau}$ of the moiré band by diagonalizing the Hamiltonian $\mathcal{H}_{\tau}$ using a plane-wave expansion. The wave function $\Phi_{\boldsymbol{k},\tau}$ contains the symmetry properties discussed above.

S2 Quantum Anomalous Hall Ferromagnets for $\Delta_{b}=\Delta_{t}$

In our study of interaction effects, we only keep the first moiré conduction band states. This approximation becomes better when the Coulomb interaction strength becomes much smaller than the band gaps to remote moiré bands. The characteristic interaction strength is set by $E_{C}=e^{2}/(\epsilon a_{M})$, where $\epsilon$ is the dielectric constant that can be controlled by the three-dimensional dielectric environment. Using $\epsilon=30$ and $\theta=1.15^{\circ}$, we find that $a_{M}\approx 12.3$ nm and $E_{C}\approx 3.9$ meV. In our moiré band structure, the band gap that separates the first and the second conduction bands is about 40 meV, which is an order of magnitude larger than the above interaction scale. The band gap that separates the first conduction band and the first valence band is the charge neutrality gap $\Delta_{\nu=0}$, which is approximately equal to $\min\{\Delta_{\bar{K}},\Delta_{\bar{K}^{\prime}}\}$. In the main text, we take $(\Delta_{b},\Delta_{t})=(30,0)$ meV, and $\Delta_{\nu=0}$ is then about 4 meV, which is, however, only slightly larger than $E_{C}$ estimated above. We note the corresponding experimental value of $\Delta_{\nu=0}$ is about 6 meV Serlin et al. .

Here we propose that $\Delta_{\nu=0}$ can be enhanced when both top and bottom hBN layers have a zero orientation angle relative to TBG. As shown in Figs. S1(d) and S2(a), $\Delta_{\nu=0}$ is enhanced to 12 meV by taking $\Delta_{b}=\Delta_{t}=30$ meV, and becomes three times the interaction scale $E_{C}(\epsilon=30)$, which makes the projection of interaction to the first conduction band a better approximation. For this TBG system with an increased $\Delta_{\nu=0}$, we perform the same analysis on interaction effects as in the main text, with results summarized in Fig. S2. We find that the main conclusions remain unchanged: (1) The quantum anomalous Hall ferromagnets (QAHF) at $\nu=3$ is generally robust against spin wave and valley wave excitations, provided that the Hartree-Fock gap is finite, as shown in Figs. S2(b) and S2(c); (2) The Curie temperature is still limited by valley wave excitations, and is reduced from the mean-field value [Fig. S2(e)]. However, there is an important difference regarding the skyrmion-antiskyrmion pair energy $\Delta_{\text{pair}}$. In particular, we find that $\Delta_{\text{pair}}$ can be lower than $\Delta_{\text{HF}}$ for certain values of $\epsilon$, as shown in Fig. S2(d). This is partially because the first moiré conduction band becomes narrower by taking $\Delta_{b}=\Delta_{t}=30$ meV. Therefore, the skyrmion-antiskyrmion pairs can be the lowest charged excitation, depending on the details of the model.

The results shown in Fig. S2 provide further theoretical evidences that the TBG QAHF are robust against particle-hole fluctuations. When the parameter $\Delta_{\nu=0}/E_{C}$ decreases, we expect that the QAHF should remain robust within a finite range of $\Delta_{\nu=0}/E_{C}$, even after effects of other bands are taken into account. An interesting question is whether there exists a critical value of $\Delta_{\nu=0}/E_{C}$, below which the QAHF becomes unstable. This question is beyond the scope of the current work, and we leave it to future study.

In Fig. S2(d), the Hartree-Fock gap $\Delta_{\text{HF}}$ at $\nu=3$ as a function of the dielectric constant $\epsilon$ has a kink around $\epsilon\approx 21$. There are similar kinks in Fig. 2(a) of the main text. The reason for this kink can be explained as follows. $\Delta_{\text{HF}}$ at $\nu=3$ is defined as $E_{+,\downarrow}^{(\text{min})}-E_{+,\uparrow}^{(\text{max})}$, where $E_{+,\downarrow}^{(\text{min})}$ is the minimum quasiparticle energy of the unoccupied spin $\downarrow$ band at valley $+K$, and $E_{+,\uparrow}^{(\text{max})}$ is the maximum quasiparticle energy of the occupied band. At small $\epsilon$ (i.e., large Coulomb interaction), $E_{+,\downarrow}^{(\text{min})}$ occurs at the center of MBZ due to a strong modification of the quasiparticle energy bands by the Coulomb interaction. By contrast, at large $\epsilon$ (i.e., small Coulomb interaction), $E_{+,\downarrow}^{(\text{min})}$ appears at the corners of MBZ. Therefore, $\Delta_{\text{HF}}$ as a function of $\epsilon$ can have a kink.