Giant orbital magnetoelectric effect and current-induced magnetization switching in twisted bilayer graphene.

Recently, quantum anomalous Hall effect with spontaneous ferromagnetism was observed in twisted bilayer graphenes (TBG) near 3/4 filling. Importantly, it was observed that an extremely small current can switch the direction of the magnetization. This offers the prospect of realizing low energy dissipation magnetic memories. However, the mechanism of the current-driven magnetization switching is poorly understood as the charge currents in graphenes are generally believed to be non-magnetic. In this work, we demonstrate that in TBG, the twisting and substrate induced symmetry breaking allow an out of plane orbital magnetization to be generated by a charge current. Moreover, the large Berry curvatures of the flat bands give the Bloch electrons large orbital magnetic moments so that a small current can generate a large orbital magnetization. We further demonstrate how the charge current can switch the magnetization of the ferromagnetic TBG near 3/4 filling as observed in the experiments.

Introduction

A bilayer graphene with a twist angle between the two graphene layers forms a quasi-two-dimensional moiré superlattice, dramatically modifying its electronic properties[1-3]. At small twist angles , the moiré potential effectively reduces the Dirac velocity[1,2] and yields flat bands at a series of magic angles[3], where electronic correlations become important. Recently, insulating[4] and superconducting phases[5], possibly driven by correlations at around 1/2 and fillings in twisted bilayer graphenes (TBG) have been observed. These discoveries have stimulated intensive theoretical and experimental works[6-36] to understand the underlying insulating and superconducting mechanisms.

More recently, experiments have unveiled the topological properties brought about by the moiré potential, as the signatures of ferromagnetism and quantum anomalous Hall effect have been experimentally observed near the 3/4 filling[37,38]. These observations are consistent with predictions[7,24] that electron–electron interactions can give rise to ferromagnetism by lifting the spin and valley degeneracy, and that quantum anomalous Hall states will be obtained when bands with a total non-zero Chern number are filled. Strikingly, these experiments have also shown that the magnetization can be switched by driving very small DC currents (from 10 to 50 nA) through the samples[37,38]. The current needed for magnetization switching is several orders of magnitude smaller than those in state of the art spin-torque devices[39]. These observations strongly suggest the possibility of realizing ultralow power magnetic memory devices in TBG. However, it is not clear how a charge current can couple to the out-of-plane magnetization of the TBG, as the charge currents in graphene layers are generally believed to be non-magnetic.

Here, we show that charge currents in TBG can induce very large orbital magnetization at general filling factors even when the sample is not ferromagnetic. We call this effect the giant orbital magnetoelectric effect. First, by symmetry analysis, we point out that due to twisting, the symmetry of bilayer graphene is reduced from D (for AB bilayer graphene) or D (for AA stacking) to D which belongs to the chiral point group. Thus, symmetry allows a magnetization to be induced by a charge current[40,41]. However, the D symmetry of TBG is still too high to allow an out-of-plane magnetization to be generated by an in-plane current for current-induced magnetization switching. Importantly, we further note that closely aligning the hexagonal boron nitride (hBN) substrate to the TBG has been essential for experimental realization of ferromagnetism and the quantum anomalous Hall state[37,38]. Including substrate-induced sublattice symmetry breaking[24-27] and strain[29,34], the symmetry of the TBG is reduced to C such that the applied current can induce a net out-of-plane magnetization[40]. Moreover, due to the large Berry curvature of the flat bands near the magic angle, the orbital magnetic moments carried by the Bloch electrons can be as large as tens of Bohr magnetons per electron even with very small strains. The lattice symmetry reduction and the large orbital magnetic moments of the electrons allow a large orbital magnetization to be induced by a small charge current. Near 3/4 filling when the Hall resistance is not quantized and the longitudinal resistance is finite (this is the experimental regime where current-induced magnetic switching has been observed[37,38]), the bulk conducting channels which carry magnetization can couple to the bulk magnetization of the sample, allowing current-controlled magnetic switching.

Results

Continuum model of strained TBG

An isolated TBG can be described by coupling the top and bottom graphene layers with a twist angle . Near the Fermi energy, the top and bottom graphene layers with the Dirac Hamiltonian at the valley can be described by a continuum model as[1-3]where denotes the Hamiltonian of the top and bottom layer respectively, is a two component creation (annihilation) operator creating (annihilating) electrons at the two A and B sublattices in the top/bottom graphene layer. The valley and the spin indices are denoted by and respectively. The momentum is defined relative to the original Brillouin zone corner that hosts the Dirac point at ,  Å is the carbon–carbon bond length[42], the rotation matrix has the form , the Fermi velocity takes the value  eV Å[42] and denotes the Pauli matrices.

An important effect of the moiré superlattice which originates from twisting is to fold the original Brillouin zone into the mini-Brillouin zones schematically shown in Fig. 1a. Both the and points of the original Brillouin zone are mapped to the mini-Brillouin zone, giving rise to fourfold degenerate minbands with both valley and spin degeneracy. In the reciprocal space, the Moiré superlattice has reciprocal vectors , , connecting the three neighboring sites of the hexagonal receiprocal lattice. The interlayer coupling is enabled when the momentum transfer between the Bloch states at different layers matches , or . The interlayer coupling Hamiltonian of the continuum model[1-3] is present in the Methods section.

Fig. 1

The energy dispersion of TBG.

a The original Brillouin zone is folded into the mini-Brillouin zone for the moiré superlattice. The solid and dashed green lines defining the large hexagons represent the original Brillouin zones of the top and bottom graphene layers respectively. b The flat bands energy dispersion from the valley with the twist angle . The dashed and solid lines are the cases with strain and , respectively. The red and blue bands represent the conduction and valence bands, respectively. The conduction flat band energy dispersion in the mini-Brillouin zone with strain for c and for d, respectively. The energy dispersion from the other valley can be mapped by the time-reversal symmetry as .

For an isolated TBG, the top and bottom graphene Hamiltonian along with the interlayer coupling respects the D symmetry[7,11,16,18,19]. At the mini-Brillouin zone corner , two massless Dirac points emerge which are protected by the composite symmetry C where is the complex conjugate operator[7,11,18,19]. However, in the two recent experiments in which a ferromagnetic state has been seen, the TBG is coupled with a hBN cladding layer aligned to the TBG to less than 1°, which empirically appears necessary to support the ferromagnetism[37,38]. In our model, the hBN substrate affects the bottom graphene layer in two aspects: (1) it breaks the C symmetry and introduces the massive gap to the Dirac Hamiltonian as shown in Fig. 1b; (2) it exerts strain on the bottom graphene layer and further reduces the crystal symmetry to C.

In this work, for simplicity we use a uniaxial strain tensor to describe the effect of strain. The strain tensor can be written aswith being the tunable parameter to characterize the strain induced displacement, the Poisson ratio for graphene[43], and the angle of the uniaxial strain relative to the zig-zag direction of the bottom graphene layer. In the presence of uniaxial strain, the real space and reciprocal space are transformed as and , respectively. Therefore the Dirac points in the bottom graphene layer are shifted from the original position to , where with being the effective gauge field from the strain[44]. By combining the sublattice symmetry breaking and uniaxial strain effect from the hBN substrate, we are able to obtain the modified bottom graphene layer Hamiltonian aswhere is the rotation matrix, the momentum is defined relative to the uniaxial strain-deformed Brillouin zone corner, and the staggered potential introduced by the hBN substrate takes  meV[45,46]. The staggered potential breaks the symmetry and reduces to symmetry, while the uniaxial strain further removes all the crystal lattice symmetry and brings the TBG down to the group.

Denoting the interlayer coupling between the strained bottom layer graphene and the unstrained top layer graphene as , the total Hamiltonian can be written aswhere is a multicomponent operator and is the Hamiltonian matrix as described in detail in the Methods section.

The energy dispersion at each band for the TBG can then be directly obtained through diagonalizing the continuum Hamiltonian in Eq. (4). For an isolated TBG, at the angle , the Hamiltonian with gives the flat bands dispersion as shown in the dashed lines of Fig. 1b. The flat bands possess two gapless Dirac points at in the mini-Brillouin zone. The energy dispersion with strain is depicted by the solid lines in Fig. 1b. The energy dispersion for the conduction band (the red band in Fig. 1b) in the whole Brillouin zone are shown in Fig. 1c and d for the unstrained and the strained cases, respectively. It is clear from Fig. 1c that the energy dispersion with spin index , valley index , and band index in general respects the C symmetry as . However, strain breaks the threefold rotational symmetry as shown in Fig. 1d, in which we set and the strain is applied along the direction of the zig-zag edge of the bottom layer graphene ( in Eq. (2)).

Orbital magnetic moment in TBG

With the Hamiltonian in Eq. (4), we can calculate the Berry curvature and the orbital magnetic moment in the out-of-plane direction of each Bloch state:and

For an isolated TBG, the Berry curvature and the magnetic moment are non-zero only at the Dirac point. In the absence of strain but in the presence of the staggered potential which is set to be  meV, the distribution of the orbital magnetic moment in the mini-Brillouin zone is shown in Fig. 2a. It respects the C symmetry as . The orbital magnetic moments are particularly large around and in the mini-Brillouin zone, where the flat band hybridizes with adjacent bands. The strength of the orbital magnetic moments can reach about with the Bohr magneton. In the presence of strain, the C symmetry of is broken as shown in Fig. 2b. If time-reversal symmetry is preserved, the orbital magnetic moment has the constraint , so that no net magnetization is allowed. However, due to the C symmetry breaking, applying a current would create an imbalance in the magnetic moment distribution and thus a net out-of-plane magnetization[40] as demonstrated in the next section.

Fig. 2

The orbital magnetic moments of the Bloch electrons.

a The orbital magnetic moments carried by the Bloch electrons in the mini-Brillouin zone with no strain. b The orbital magnetic moments of the electrons when a uniaxial strain characterized by is introduced. The staggered potential is set to be  meV in both cases. In b, C symmetry is broken and the Brillouin zone is deformed.

Magnetoelectric response in TBG

In quasi-two-dimensional materials with finite magnetoelectric response, the electric field induced magnetization can be described aswith , and the magnetoelectric susceptibility. As shown in refs. [40,41], the general forms of the components of the magnetoelectric susceptibility tensor can be determined by the crystal symmetry of the material. The general forms of for point groups D, C, and C which are relevant to TBG are shown in Table 1. It is clear from Table 1 that it is possible to generate an out-of-plane magnetization by in-plane electric fields only if the crystal point group symmetry is reduced to C.

Table 1

Magnetoelectric susceptibility pseudotensor for D, C, and C point group.

Point group	α	Point group	α
D_6	α_∥000α_∥000α_zz	C_3	α_∥−α^−0α^−α_∥000α_zz
C_1	α_xxα_xyα_xzα_yxα_yyα_yzα_zxα_zyα_zz

with are in general the elements in . In D and C, is denoted as . In C, the antisymmetric off diagonal element is denoted as .

To calculate for TBG, we can use the linear response theory which gives[47,48]where , is the Fermi Dirac distribution function, is the group velocity, is the effective scattering time, and the total magnetic moment is composed of both the orbital magnetic moment and the spin magnetic moment with the Lande g factor .

To be specific, we apply a uniaxial strain with along the zig-zag edge direction of the bottom layer graphene. The orbital magnetization in the Brillouin zone in the presence of strain is shown in Fig. 2b. The resultant magnetoelectric susceptibility can then be evaluated assuming the electron scattering time to be  ps[49]. For the conduction band , the magnetoelectric susceptibility is shown in Fig. 3a as a function of the Fermi energy, where the Cartesian coordinate is set to have the -axis along the angular bisector between the two zig-zag directions of the top and bottom graphene layers. The magnetoelectric susceptibility is maximized near the energy with the largest density of states. Interestingly, are still very large even when the density of states is very low. This is because the orbital magnetizations carried by the Bloch states near are very large as a result of the Berry curvatures of the flat bands. This allows a large magnetization to be induced by a small current. As shown in the Supplementary Fig. 2, the current-induced orbital magnetization can be even stronger when strain is increased.

Fig. 3

The magnetoelectric response in strained TBG.

a The magnetoelectric susceptibilities , , and the density of states (DOS), both as a function of the Fermi energy from bottom to top of the conduction band. We have set . b The induced magnetization at the electric field strength  V/m along different in-plane angles. Increasing radius in the polar plot denotes the Fermi energy increasing from the conduction band bottom to the top. The twist angle is set to be .

Assuming an external electric field of  V/m, we obtain the out-of-plane magnetization under different electric field directions as shown in Fig. 3b, where the increasing radius in the polar plot denotes the Fermi energy increases from the conduction band bottom to the top. The magnetization can reach , 1–2 orders larger than in the largest Rashba spin–orbit coupling materials such as Au (111) surfaces and Bi/Ag bilayers[50,51]. The current-induced magnetization is anisotropic with respect to the direction of the current and it switches sign under reversal of the electric field. It is important to note that the current-induced magnetization discussed here can appear at a general filling factor even absent spontaneous ferromagnetism in the sample. This current-induced magnetization should be observable experimentally through optical Kerr effects as in the case of transition metal dichalcogenides[52].

Current-induced magnetization switching in TBG

TBG in the non-interacting limit possess valley and spin degeneracy for each flat band[1-3,7,11,18,19]. However, near the magic angles, the narrow band width at the Fermi level magnifies the role of interactions, and interaction-driven spontaneous symmetry breaking is observed experimentally[37,38]. Specifically, at 3/4 filling of the conduction band in hBN-aligned TBG with inter-graphene twist angle (ref. [37]) a giant anomalous Hall effect of order has been reported; and for TBG with twist angle [38], quantized anomalous Hall effect has been reported, in both cases at zero external magnetic field. Hysteresis in the Hall conductance under out-of-plane magnetic fields suggests spontaneous ferromagnetism with out-of-plane magnetization.

The presence of net magnetization as revealed by anomalous Hall resistance[37,38] indicates that the spin and/or valley degeneracies are lifted, possibly by interactions[7,24-27]. As a result, there are four bands (which originated from the fourfold degenerate conduction band in Fig. 1b) labeled by the spin indices and valley indices available for the electrons to fill as shown in Fig. 4a and b. To take into account the simplest possible spin and valley polarization phenomenologically, the dispersion of the four bands is written as where describes the original conduction band dispersion from valley without interaction and is the spin- and valley-dependent energy shift due to interactions.

Fig. 4

Current induced magnetization switching.

The interaction-renormalized bands at 3/4 filling for heterostrain are shown in a for fully-filled spin- and valley-polarized bands, and in b for spin-polarized but valley-unpolarized bands. c Free energy as a function of magnetization for several values of applied electric field along the direction. At the coercive electric field , one local minimum of the free energy collapses and only one minimum remains. d The magnetic hysteresis curve induced by the electric field along the direction. The coercive field is estimated to be around  V/m.

At filling factor 3/4, if the three bands with lower energy are completely filled as depicted in Fig. 4a, the TBG should display the quantum anomalous Hall effect. At the same filling factor 3/4, the top two bands could instead each be partially filled as seen in Fig. 4b, in which case the TBG would have a bulk conducting channel in parallel with the anomalous Hall conductance. The scenario of Fig. 4b may be a good representation of experiments where the Hall conductance is not quantized and bulk conducting channels exist[37].

To connect our theory with experiments, we note that the spontaneous ferromagnetism in TBG can be described by the Landau’s free energy density asBelow the critical temperature, , , generating a finite magnetization order parameter at and the magnetic susceptibility reads . In the presence of external magnetic field , the magnetization switches sign at the coercive magnetic field . Note that and can be obtained once and are calculated using the continuum model introduced previously with the energy of the bands shifted by . Given , the total magnetization and the magnetic susceptibility can be evaluated aswhere is the z-component of the total magnetic moment of a Bloch wavefunction of the flat bands. In the partially polarized state shown in Fig. 4b with  meV, we find that , and the coercive magnetic field  mT.

To understand the coupling between the electric field and the magnetic field, we note that the total magnetization is changed to where is the magnetization induced by the current. As a result, the Landau free energy in the presence of an electric field can be written aswhich clearly shows that the magnetization of the sample couples to the electric field. Figure 4c depicts the free energy landscape as a function of magnetization changes for different electric field strength, using realistic parameters. By assuming the current is passed in the y-direction and by calculating , the resulting hysteresis loop of magnetization as a function of electric field is determined. The minimal electric field needed to switch the magnetization is estimated to be about 113 V/m. In a recent experiment[37], the longitudinal resistance is measured to be and the length between the contacting leads is estimated to be 5 μm. As a result, the coercive electric field at  V/m gives the coercive DC current  nA, which matches well with the experimental values of  nA[37]. Since many of the details such as the strain, the band structure of the sample, the shifts of the polarized bands, etc. will affect the coercive current, the specific value of the coercive electric field calculated here can only be a rough estimation.

Discussion

In the above sections, using a continuum model of TBG and incorporating the effects of sublattice symmetry breaking and strain, the magnetoelectric response was calculated. Here, we would like to emphasize that the analysis based on symmetry is very general. The exact form of the strain is not important. The breaking of the D symmetry can come from other sources such as spatial inhomogeneity in the chemical potential or twist angles. The detailed source of symmetry breaking will not affect our conclusion that currents can induce magnetization in TBG. Moreover, the current-induced magnetization effect can appear even when the system itself is not ferromagnetic (for example, in the absence of valley polarization). Therefore, we expect that other materials with low crystal symmetries such as twisted bilayer-bilayer graphene[53-55], twisted hBN-graphene heterostructure[56,57], twisted transition metal dichacolgenides[58], and gapped bilayer graphene[59] with strain will exhibit similar magnetoelectric effects, although the magnitude of the magnetoelectric response will depend on the details of the materials. The current-induced orbital magnetization predicted can be tested by magneto-optical Kerr effect in experiments[52].

Another important point is that in the experimental regime where current-induced magnetization switching is demonstrated, the Hall resistance R is not quantized and the longitudinal resistance R is finite[37,38]. The currents can flow between domains with different magnetization. As the symmetry of the problem is still C even including the domains, the bulk currents can carry out-of-plane magnetization and switch the magnetization of the domains. However, a calculation incorporating domains is beyond the scope of the current study.

Our picture of current-induced magnetic switching does not apply directly to quantum anomalous Hall states with an insulating bulk when the current is carried by the edge states only. To obtain the current-induced magnetization, we assumed that the scattering time () in the system is finite as shown in Eq. (8). This assumption does not apply to chiral edge states. Serlin et al.[38] argued that even edge states which do not carry net magnetization can also switch the direction of the magnetic domains, giving an effect proportional to where is the current carried by the edge states. In contrast, in the present work, the magnetoelectric effect of the bulk currents couples the electric field linearly to the magnetization as shown in Eq.  (12).

It is also worth noting that the current-induced magnetization in TBG is purely orbital in nature. It is different from the magnetoelectric effect induced by spin–orbit coupling in noncentrosymmetric materials[60,61] studied previously. It is also interesting to note that the orbital magnetization can be strongly affected by strain. In this work, we only discussed the strain induced naturally by the hBN substrate. Experimentally, one can induce a much larger strain on the TBG artificially. In this case, the current-induced magnetization could be further enhanced. The orbital magnetization of some of the Bloch states in the Brillouin zone can even reach a hundred Bohr magnetons with moderate strain as shown in the Supplementary Fig. 1. In this case, even larger orbital magnetoelectric effects could be realized in TBG.

Methods

Interlayer coupling Hamiltonian for the TBG

In the continuum model description[1-3], the state at from one layer will couple with the state at from the other layer if matches , , or , so the interlayer coupling Hamiltonian readswith the tunneling matrixIn the presence of uniaxial strain in the bottom layer graphene, the reciprocal vectors are deformed as , , and the tunneling matrix are modified as , , , where the detailed forms are presented in the Supplementary Note 1.

The Hamiltonian matrix for TBG coupled with hBN substrate

The Hamiltonian for the TBG on a hBN substrate readswhere has infinite components representing the series of states , with . For example, in the truncated basis of , the Hamiltonian matrix has the form:withIn our calculation, we consider 42 sites in the hexagonal reciprocal lattice and is a matrix in the truncated basis. The Hamiltonian matrix for the TBG respects the time-reversal symmetry aswith the energy dispersion for the band index , valley index , spin index .